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Abstract 


New  families  of  highly  non-Keplerian  orbits  can  be  generated  using  continuous  low-thrust  propulsion. 
These  include  orbits  displaced  above/below  the  orbit  plane  of  a  conventional  Keplerian  orbit  using  out-of- 
plane  thrust,  or  inside/outside  a  Keplerian  orbit  using  radial  thrust.  For  a  given  non-Keplerian  orbit  radius, 
displacement  distance,  and  orbit  period  the  required  thrust  magnitude  and  direction  can  be  determined  an¬ 
alytically,  along  with  the  orbit  stability  properties. 

A  key  open  issue  in  the  analysis  of  the  families  of  non-Keplerian  orbits  is  the  generation  of  a  mapping 
from  the  non-Keplerian  orbit  radius,  displacement  distance  and  orbit  period  to  a  corresponding  set  of  oscu¬ 
lating  Keplerian  orbital  elements.  This  is  of  importance  both  to  support  the  future  operational  use  of  non- 
Keplerian  orbits  and  to  understand  the  signature  of  a  non-Keplerian  orbit  as  seen  in  spacecraft  tracking  data. 
Since  families  of  non-Keplerian  orbits  are  generated  using  a  strong,  continuous  perturbation,  their  osculat¬ 
ing  Keplerian  orbital  elements  will  be  time  varying,  although  still  periodic.  Similarly,  the  inverse  problem 
is  also  of  importance  to  deduce  the  non-Keplerian  orbit  radius,  displacement  distance,  and  orbit  period  from 
observations  of  Keplerian  elements.  Once  these  elements  are  known,  the  thrust-induced  acceleration  mag¬ 
nitude  and  direction  required  for  the  non-Keplerian  orbit  can  also  be  determined. 

In  order  to  pursue  these  questions,  the  key  objectives  of  this  work  are,  therefore,  to:  (1)  map  the  proper¬ 
ties  of  families  of  highly  non-Keplerian  orbits  onto  the  classical  Keplerian  elements,  (2)  generate  an  inverse 
mapping  from  Keplerian  elements  to  the  properties  of  the  non-Keplerian  orbit,  and  (3)  determine  the  key 
signatures  of  non-Keplerian  orbits  in  Keplerian  element  tracking  data.  These  objectives  represent  a  largely 
unexplored  aspect  of  the  dynamics  of  families  of  non-Keplerian  orbits  and  offer  a  route  towards  their  future 
operational  use. 

This  report  presents  a  mapping  between  the  elements  of  highly  non-Keplerian  orbits  and  classical  orbital 
elements.  Three  sets  of  elements  are  discussed  and  mappings  are  derived  in  closed,  analytical  form  for  both 
the  direct  and  inverse  problem.  Advantages  and  drawbacks  of  the  use  of  each  set  of  elements  are  discussed. 
A  sensitivity  analysis  is  performed  on  all  mappings  presented.  The  spacecraft  thrust-induced  acceleration 
used  to  generate  families  of  non-Keplerian  orbits  is  extracted  from  the  inverse  mapping  from  the  osculating 
orbital  elements.  The  key  signatures  of  highly  non-Keplerian  orbits  in  Keplerian  element  tracking  data  are 
determined  through  a  set  of  representative  test  cases. 
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1.  Highly  Non-Keplerian  Orbits:  Introduction 

In  classical  Keplerian  mechanics,  an  object  of  negligible  mass  orbits  around  an  attractor  due  to  the  gravitational 
interaction  existing  between  the  two  bodies.  Celestial  mechanics  traditionally  deals  with  Keplerian  orbits,  treating  any 
non-Keplerian  effect  (e.g.  solar  radiation  pressure,  atmospheric  drag,  gravitational  asymmetries  due  to  the  Earth’s 
oblateness,  third-body  interaction)  as  a  perturbation  to  the  classical  Keplerian  orbit.  However,  a  non-Keplerian  orbit 
is  defined  as  an  orbit  in  which  a  perturbative  or  propulsive  acceleration  acts  in  addition  to  the  gravitational  attraction 
of  the  primary  body.  Following  the  definition  given  in  McKay  et  al.,1  a  highly  non-Keplerian  orbit  is  characterized  by 
the  magnitude  of  the  average  propulsive  acceleration  over  one  orbit  || aav  ||  equal  or  greater  than  the  sum  of  gravitational 

and  centripetal  force  experienced  by  the  orbiting  body  ||V£/|| .  That  is,  introducing  a  parameter  A  to  represent  this 

ratio,  a  highly  non-Keplerian  orbit  is  characterized  by  A  =  ||«av||/||Vf/||  ~  1 .  For  the  sake  of  conciseness,  and  because 

only  Keplerian  and  highly  non-Keplerian  orbits  are  considered  in  this  paper,  all  the  highly  non-Keplerian  orbits  will 
be  simply  referred  to  as  non-Keplerian  orbits  (NKOs)  in  this  report. 

New  families  of  NKOs  can  be  generated  using  continuous  low-thrust  propulsion.1'4  These  include  orbits  displaced 
above/below  the  orbit  plane  of  a  conventional  Keplerian  orbit  using  out-of-plane  thrust,  or  inside/outside  a  Keplerian 
orbit  using  radial  thrust.  For  a  given  NKO  radius,  displacement  distance  and  orbit  period  the  required  thrust  magnitude 
and  direction  can  be  determined  analytically.2  In  addition,  the  linear  stability  properties  of  the  families  of  non-Kep¬ 
lerian  orbits  can  be  determined  and  it  can  be  shown  that  while  unstable  families  of  orbits  exist  (saddle  x  center  eigen¬ 
value  spectrum)  they  are  in  principle  controllable  (full  rank  controllability  matrix)  using  feedback  to  the  thrust  direc¬ 
tion  and/or  modulation  of  the  thrust  magnitude.3  Moreover,  the  nonlinear  stability  properties  of  NKOs  have  also  been 
investigated,  and  the  sufficient  and  necessary  conditions  for  stability  of  the  motion  near  the  equilibria  have  been  ob¬ 
tained.5 

In  addition,  by  terminating  the  low-thrust  perturbation,  a  Keplerian  orbit  will  result  allowing  Keplerian  and  non- 
Keplerian  orbits  to  be  patched  together  to  generate  rich  and  complex  families  of  composite  orbits  that  have  yet  to  be 
fully  explored.6  Potential  applications  of  NKOs  include  orbits  displaced  above/below  the  geostationary  ring  to  increase 
the  number  of  available  slots  for  communications  platforms4  and  on-orbit  inspection  by  formation-flying  above/below 
or  inside/outside  the  orbit  of  a  target  spacecraft.7, 8  For  example,  a  displacement  distance  of  35  km  north/south  of  the 
geostationary  ring  allows  a  platform  on  a  displaced  highly  non-Keplerian  orbit  to  sit  above/below  the  standard  70  km 
station-keeping  box  of  a  conventional  geostationary  platform.  Such  an  orbit  would  require  a  thrust  of  200  mN  for  a 
1000  kg  spacecraft,  which  is  achievable  with  a  single  QinetiQ  T6+  thruster.4 

To  date,  several  works  have  been  carried  out  to  study  non-Keplerian  orbits.  As  well  shown  in  the  survey  presented 
by  McKay  et  al.,1  the  problem  has  been  considered  from  several  aspects  by  taking  into  account  different  propulsion 
technologies  (e.g.  impulsive  thrust,6  solar  electric  propulsion,9  solar  sail,10, 11  and  hybrid  SEP/sail4, 12)  and  different 
dynamical  models  (e.g.  two-body,2, 3  circular  restricted  three -body, 13  elliptical  restricted  three-body,14  and  restricted 
four-body15).  Different  mission  applications  are  also  outlined  for  non-Keplerian  orbits,  such  as  telecommunication,7 
Earth  observation,16  planetary  science,17  and  geoengineering.18  More  recently,  Wang  et  al.19, 20  analyzed  the  relative 
motion  between  heliocentric  displaced  NKOs.  Describing  the  motion  by  means  of  classical  Keplerian  elements  and 
modified  equinoctial  elements,21  they  provided  an  analytical  solution  to  the  relative-motion  problem. 

As  discussed  above,  all  prior  studies  on  non-Keplerian  orbits  have  focused  on  orbit  and  mission-design.  A  key 
open  issue  in  the  analysis  of  the  families  of  non-Keplerian  orbits  is  the  generation  of  a  mapping  from  the  NKO  geom¬ 
etry  to  a  corresponding  set  of  osculating  orbital  elements.  This  is  of  importance  both  to  support  the  future  operational 
use  of  non-Keplerian  orbits  and  to  understand  the  signature  of  a  non-Keplerian  orbit  as  seen  in  spacecraft  tracking 
data.  Since  families  of  NKOs  are  generated  using  a  strong,  continuous  perturbation,  their  osculating  orbital  elements 
will  be  time  varying,  although  still  periodic.  Similarly,  the  inverse  problem  is  also  of  importance  to  deduce  the  NKO 
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geometry  from  observations  of  orbital  elements.  Once  these  elements  are  known,  the  thrust-induced  acceleration  mag¬ 
nitude  and  direction  required  for  the  non-Keplerian  orbit  can  also  be  determined. 

In  order  to  pursue  these  questions,  the  key  objectives  of  this  study  are,  therefore,  to: 

7.  Map  the  properties  of  families  of  highly  non-Keplerian  orbits  onto  the  classical  orbital  elements 

2.  Generate  an  inverse  mapping  from  orbital  elements  to  the  properties  of  the  non-Keplerian  orbit 

3.  Determine  the  key  signatures  of  non-Keplerian  orbits  in  orbital  element  tracking  data 

These  objectives  represent  a  largely  unexplored  aspect  of  the  dynamics  of  families  of  non-Keplerian  orbits  and 
offer  a  route  towards  their  future  operational  use. 

The  report  is  divided  into  two  main  Sections.  Section  2  describes  the  case  of  a  simplified  model,  in  which  the  NKO 
orbital  plane  is  always  parallel  to  the  equatorial  plane.  On  the  other  hand,  a  generalized  model  with  no  assumptions 
on  the  NKO  orbital  plane  is  described  in  Section  3.  In  the  case  of  the  simplified  model,  the  forward  mappings  from 
NKO  geometry  to  osculating  orbital  elements  are  discussed  in  Section  2.1.  The  description  of  the  mappings  and  their 
sensitivity  analyses  are  analytically  investigated  in  Sections  2.1.1  and  2.1.2,  respectively,  whereas  Section  2.1.3  shows 
the  numerical  test  cases  taken  into  account  for  the  study  of  the  mappings.  In  Section  2.2,  the  inverse  mappings  from 
osculating  orbital  elements  to  non-Keplerian  orbits  are  discussed  together  with  their  sensitivity  analyses.  The  thrust- 
induced  acceleration  needed  to  maintain  the  chosen  non-Keplerian  orbits  is  discussed  in  Section  2.3.  In  the  case  of  the 
generalized  model,  the  analytical  description  of  the  forward  mappings  is  investigated  in  Section  3.1.  Sections  3.1.1 
and  3.1.2  provide  the  analytical  derivation  of  the  forward  maps  by  using  the  classical  Keplerian  elements  and  the 
modified  equinoctial  elements,  respectively.  Numerical  test  cases  are  shown  and  discussed  in  Section  3.1.3,  whereas 
a  summary  of  the  signatures  of  a  spacecraft  being  forced  on  a  NKO  is  provided  in  Section  3.1.4.  A  study  on  the 
generalized  inverse  mapping  is  given  in  Section  3.2  which  provides  an  analytical  description  of  all  families  of  NKOs 
that  share  the  same  Cartesian  state.  Lastly,  the  impact  of  noise  is  discussed  in  Section  3.3.  Conclusions  and  final 
remarks  are  drawn  in  Section  4. 
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2.  Simplified  Model:  Orbital  Plane  Parallel  to  Equatorial  Plane 

Non-Keplerian  orbits  can  be  obtained  by  considering  the  dynamics  of  a  low-thrust  propelled  spacecraft  in  a  rotating 
reference  frame.2  The  angular  velocity  of  rotation  of  the  reference  frame  m  ,  together  with  the  out-of-plane  displace¬ 
ment  z  and  the  radius  of  the  orbit  p  ,  are  used  as  free  parameters  of  the  problem.  Stationary  solutions  of  the  equations 
of  motion  in  this  rotating  reference  frame  correspond  to  periodic,  displaced,  circular  orbits  with  the  orbital  plane 
parallel  to  the  X—Y  plane  when  viewed  from  an  inertial  reference  frame.  In  this  paper,  the  Earth-Centered  Inertial 
frame  jx,F,zj  is  considered.  The  acceleration  needed  to  generate  such  stationary  solutions  can  be  derived  in  a 

closed,  analytical  form.2  The  required  thrust  vector  lies  in  the  plane  spanned  by  the  radius  vector  and  the  axis  of 
symmetry,  as  shown  in  Figure  1 .  Therefore,  the  propulsive  acceleration  vector  can  be  described,  as  a  function  of  the 

NKO  elements  xNKO  -  [z,p,zuf  ,  by  its  magnitude  ||a||  and  pitch  angle  a  ,  as  follows. 

\\a(z,p,v)\\  =  ^7 

< 

tana(z,p,nj)  =  — 

z 

The  term  m*  =  ^p/r3  in  Eq.  (2.1)  is  the  orbital  angular  velocity  of  a  Keplerian  orbit  of  radius  r  =  *Jp2  +z2  . 

Viewed  from  an  inertial  reference  frame,  the  orbits  generated  by  the  acceleration  described  in  Eq.  (2.1)  correspond  to 
circular  orbits  displaced  above  the  central  body.  The  angular  velocity  of  the  rotation  of  the  reference  frame  zu  corre¬ 
sponds  to  the  angular  velocity  of  the  circular  orbit  viewed  from  an  inertial  reference  frame. 


Y 


+z2mi 


{  \2 

m 


(2.1) 


Figure  1.  Displaced  orbit  with  thrust-induced  acceleration  for  the  simplified  model. 

The  acceleration  defined  in  Eq.  (2.1)  is  that  required  to  generate  stationary  solutions  in  the  rotating  reference  frame 
(i.e.  ||a ||  =  \\VU\\ ).  Therefore,  such  solutions  are  characterized  by  ?i  =  1  and  so  are  defined  as  highly  non-Keplerian 
orbits. 

Three  families  of  NKOs  can  be  defined  based  on  the  value  of  the  angular  velocity  of  the  rotating  reference  frame.2 
Type  1  NKOs  are  defined  as  those  orbits  with  the  minimum  required  acceleration.  From  Eq.  (2.1),  the  requirement  of 
minimum  acceleration  leads  to  an  angular  velocity  of  the  rotating  reference  frame  for  Type  1  NKOs  defined  as  shown 
in  Eq.  (2.2).  A  second  family  of  NKOs  is  characterized  by  orbits  synchronous  with  a  body  on  a  circular  Keplerian 
orbit  in  the  z  =  0  plane  with  the  same  orbit  radius.  Lastly,  a  third  family  of  NKOs  is  defined  by  setting  the  orbital 
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period  to  a  fixed  value.  The  angular  velocities  that  characterize  the  aforementioned  three  families  of  NKOs  are  shown 
in  Eq.  (2.2). 

®Typel  =4^’  ®Type2  =  ^//?\  ®'TWe3=®'o  (2-2) 

In  this  Section,  the  mappings  between  NKO  elements  and  the  osculating  orbit  are  derived  in  closed,  analytical 
form  for  both  the  direct  and  inverse  problem.  Therefore,  the  formulation  of  the  spacecraft  thrust -induced  acceleration 
given  in  Eq.  (2.1)  is  derived  as  a  function  of  the  osculating  orbital  elements. 


2.1.  Forward  Map:  Non-Keplerian  Orbit  to  Osculating  Orbital  Elements 

Three  sets  of  elements  are  chosen  to  describe  the  osculating  Keplerian  orbits: 

xkep  ~  [a> e >  U  Q  <$■>  z9]r  Classical  Keplerian  elements  (KEP) 

x MEE  =  \p,f,g,h,k9 Lf  Modified  equinoctial  elements  (MEE)  (2.3) 

xaiom  -  [ hT ,  eT ,  Lj  Augmented  integrals  of  motion  (AIoM) 

From  Mclnnes,2  the  three  families  of  NKO  are  fully  described  by  the  following  three  elements  in  an  Earth-centered 
cylindrical  reference  frame. 


p  Radius  of  the  orbit 
z  Out-of-plane  displacement 
uj  Angular  velocity  of  the  spacecraft 

Therefore,  the  Cartesian  position  vector  r  is  simply  given  by 

pcosipnt) 

z 


(2.4) 


(2.5) 


and  the  Cartesian  velocity  vector  v  is 


-67/9  sin  (W) 
mp  cos  [mt^ 
0 


(2.6) 


Because  the  NKOs  considered  are  circular,  it  is  worth  noting  that  the  radial  component  of  the  velocity  is  always  zero, 
as  it  can  be  easily  verified  by  computing  vr  =  v  •  r  =  0 . 


The  orbital  angular  momentum  h  can  be  computed  by  using  the  aforementioned  expressions  of  position  and  ve¬ 
locity  as 


h  -  r  xv 


—jupz  cos  (rut ) 
—zupz  sin  {jut ) 
mp2 


(2.7) 


The  analytical  formulation  of  the  three  mappings  is  detailed  in  Section  2.1.1,  whereas  the  sensitivity  analysis  per¬ 
formed  on  the  mappings  is  described  in  Section  2.1.2.  Section  2.1.3  provides  numerical  test  cases.  Lastly,  the  three 
mappings  are  compared  through  their  advantages  and  drawbacks  in  Section  2.1.4. 
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2.1.1.  Mappings 

In  the  following  sub-sections,  the  forward  mappings  from  non-Keplerian  orbit  to  osculating  orbital  elements  are 
derived  for  the  three  sets  of  elements  considered. 


2.I.I.I.  Classical  Keplerian  Elements 


Starting  from  the  Cartesian  state  vector  (Eqs.  (2.5)  -  (2.6)),  the  mapping  from  NKO  elements  \z,p,tu ]r  to  oscu- 
lating  classical  Keplerian  elements  is  obtained,  as  described  in  Algorithm  4.2  from  Reference  22. 

1.  Inclination  ( i ) 


i  =  arccos 


(AZ)  = 


arccos 


4777 


(2.8) 


2.  Right  Ascension  of  the  Ascending  Node  ( Q  ) 

The  line  of  the  nodes  shall  be  defined  in  order  to  find  the  expression  for  the  right  ascension  of  the  ascending  node 
(RAAN).  The  line  of  the  nodes  is  defined  as 


jv=iiM 

Z  xh 


sin(sj/L)sgn(z) 

-cos(^t)sgn(z) 

0 


(2.9) 


It  is  important  to  note  that  the  line  of  the  nodes  is  not  defined  in  case  of  an  in-plane  orbit  (i.e.  z  =  0 )  because 
sgn(O)  =  0  .  Therefore,  the  line  of  the  nodes,  in  this  case,  is  arbitrarily  chosen  to  point  in  the  X  direction.  Therefore, 
the  general  formulation  of  the  line  of  the  node  is 


N  = 


cos(Q) 

sin(Q) 

0 


1 

0 

0 

sin(^7/L)sgn(z) 

-cos(670sgn(z) 


0 


if  z  =  0  a  jn  -  m2 p3 


otherwise 


(2.10) 


Defining  the  four-quadrant  inverse  tangent  as*  arctan (X , Y )  =  tan  1  (Y/X)  +  (7r/ 2) sgn (F)  (l  - sgn (X ))  ,  the  right  as¬ 
cension  of  the  ascending  node  is  simply  given  by  Eq.  (2.11),  in  which  Q*  =  arccos  (sin  (mt)  sgn  (z))  . 


Q  =  arctan •  X,N  -  Y^  = 


0 

Q* 

27T-C1* 


if  z  =  0 

if  (cOs(ftTt)  <  0  A  z  >  o)  v(cos(ft7t)  >  0  A  z  <  o) 
otherwise 


(2.11) 


*  Definition  from  MATLAB  User’s  Guide,  Ver.  R2014b,  atan2  -  Symbolic  four- quadrant  inverse  tangent. 
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3.  Eccentricity  (  e ) 

The  eccentricity  vector  is  given  by 


A  v  xh 

e  =  —r  H - 


(2.12) 


Using  the  vector  triple  product  rule  and  remembering  that  vr  =  0 ,  the  expression  for  the  eccentricity  vector  in  Eq. 
(2.12)  can  be  rewritten  as 


e  = 


2 

m  p 


p  4777 


pcos(mt) 

psin(mt) 


(2.13) 


The  value  of  the  eccentricity  is,  therefore,  given  by  the  norm  of  the  eccentricity  vector,  as 


__2  2 
-  UJ  p 


l 


p1  +z2 


(2.14) 


Note  that  a  circular  planar  Keplerian  orbit  is  characterized  by  p  =  m1  p3 ,  which  satisfies  the  zero-eccentricity 
condition.  However,  the  eccentricity  vector  is  not  defined  in  this  case  and  this  is  an  issue  in  the  following  computation 
of  the  true  anomaly.  Therefore,  the  eccentricity  unit  vector,  in  this  case,  is  arbitrarily  chosen  to  point  in  the  X  direc¬ 
tion.  Therefore,  the  general  formulation  of  the  eccentricity  vector  is 


e 


1 

0 

0 


_2  2 

UJ  p 


\ 

p  cos  ^zut  ^ 

1 

//  o  o 

psmfnjt) 

4p  +z  > 

z 

if  z  =  0  a//  =  m2p3 


otherwise 


(2.15) 


4.  Semi-major  axis  (a) 

The  semi -major  axis  is  given  by 


H2 

m\Ip2+z2 

//(i-e2) 

2  ju-m2p2s]p2  +z2 

(2.16) 


5.  Argument  of  perigee  ( co ) 

The  argument  of  perigee  is  given  by 

co  =  arccos^iV  -e^j  (2.17) 
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After  some  mathematical  manipulations  and  arbitrarily  choosing  co  =  0  in  the  case  of  planar  circular  Keplerian 
orbit,  the  general  expression  for  the  argument  of  perigee  is 


co-< 


0 

if  z  =  0  a//  =  m2p3 

mt 

if  z  =  0  a//  <  m2p3 

n  +  mt 

if  z  =  0  a//  >  m2p3 

n /  2 

if  |z  >  0  a//  <  m2 p‘ 

3^/2 

if  ( z  >  0  a//  >  m2 p‘ 

(2.18) 


6.  True  anomaly  (  3 ) 

The  true  anomaly  is  given  by 


cos 


(*)= 


e  r 


(2.19) 


In  the  case  of  a  circular  planar  Keplerian  orbit,  Eq.  (2.19)  becomes  Eq.  (2.20),  because  of  the  choice  of  the  eccen¬ 
tricity  vector  shown  in  Eq.  (2.15). 


cos 


(i9)  =  cos  ( mt ) 


3-mt 


(2.20) 


On  the  other  hand,  a  NKO  is  characterized  by  the  expression  of  the  true  anomaly  given,  after  some  mathematical 
manipulations,  as 


3  - 


[o  if  fi  <  m2p2^Jp2  +  z 


\n  if  p>m2  p2  p 2  +  z 

Therefore,  the  general  expression  for  the  true  anomaly  is  found  by  unifying  Eqs.  (2.20)-(2.21),  as  follows. 


(2.21) 


3  =  { 


mt  if  z  =  0  A/j  =  in1  p2' 

0  if  ju  <  m2p2^jp2  +z2 


(2.22) 


n  if  ju  >  m2p2^p2  +z2 

Note  that  the  case  ^  0  a/j  =  m2p2*Jp2  +z2  j  is  not  taken  into  account  inside  Eq.  (2.22).  However,  this  case 

corresponds  to  an  osculating  circular  orbit  and  this  is  not  possible  for  the  types  of  NKOs  taken  into  account  within 
this  study. 
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Summary 

The  forward  map  from  NKOs  to  osculating  classical  Keplerian  elements  described  in  Section  2. 1.1.1  is  summa¬ 
rized  as: 


mIp2  +  z2 


2//-E7  2  p2  yfp1  +  Z' 

ju-sr2p2ijp2  +z2 


i  -  arccos 


Q  =  < 


arccos 


(sin(W)sgn(z)) 


if  z  -  0 

if  (cos  [mt )  <  0  a  z  >  0)  v  (cos  [mt)  >  0  a  z  <  0) 


co  =  < 


3  = 


2 n  -  arccos  (sin  [mt)sgn  (z))  otherwise 

0  if  z  =  o  A  jU  =  m2  p2 

mt  if  z  =  0  a  ju  <  m2  p2 

7t  +  mt  if  z  =  0  a  jU  >  m2  p2 

nj 2  if  >  0  aju  <  m2  p2^p2  +z2  j  v|z  <  0  a//  >  m2p2^p2  +z2  j 

3tt/2  if  >  0 /\p  >  m2p2^p2  +  z2  jv|z  <  0 a p<  m2p2^Jp2  +  z2  j 

mt  if  z  =  0  a//  =  m2p2 
0  if  p  <  m2p2*Jp2  +z2 
n  if  p  >  m2 p2  ^ p 2  +z2 


(2.23) 
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2.I.I.2.  Modified  Equinoctial  Elements 

In  order  to  have  a  mapping  with  no  singularities,  a  map  from  Cartesian  position  and  velocity  to  Modified  Equi¬ 
noctial  Elements  (MEE)  is  derived.  MEE  have  been  introduced  to  avoid  singularities  occurring  in  Keplerian  orbits  in 
the  case  of  planar  and/or  circular  orbits  and  are  defined  as  follows:21 

p  =  a(l-e2^j 
f  -  ecos(&>  +  Q) 

g  =  esin(®  +  Q)  (2.24) 

h  =  tan  (//2)  cos  (Q) 
k  =  tan(//2)sin(Q) 

L  —  co  +  £2  + 19 

Starting  from  the  expressions  of  the  Cartesian  position  and  velocity  vectors  (Eqs.  (2.5)  -  (2.6))  and  the  analytical 
mapping  of  the  osculating  classical  Keplerian  elements  described  in  Section  2. 1.1.1  (Eq.  (2.23)),  the  mapping  from 

NKO  elements  [z,p,cuf  to  osculating  modified  equinoctial  elements  is  derived  in  the  following  sub-sections. 

1.  Semi-latus  rectum  (  p ) 


||ftf  rp!(p’+z’) 
p  p 


(2.25) 


2.  In-plane  element  (  f  ) 


In  order  to  derive  the  general  formulation  for  the  in-plane  element  /  ,  let  us  first  consider  the  case  of  zero-dis¬ 
placement  orbits  (i.e.  z  =  0 ). 

From  Eq.  (2.1 1),  the  value  of  the  right  ascension  of  the  ascending  node  for  a  zero -displacement  orbit  is  simply 

Q|  _o  =  0  (2.26) 


From  Eq.  (2.18),  the  value  of  the  argument  of  perigee  for  a  zero -displacement  orbit  is 


0 


co-{ 


cut 


K  +  CUt 


if  z  =  o  Ap  =  cu2p3 
if  z  =  0  /\p  <  cu2 p3 
if  z  =  0  Ap  >  cu2 p3 


From  Eqs.  (2.26)-(2.27),  the  value  of  the  longitude  of  perigee  cd  +  Q.  is 


0 


co  +  —  < 


cut 


7C  +  CUt 


if  z  =  0  Ap  =  cu2  p3 
if  z  =  0  a  p  <  cu2 p3  => 

if  z  =  0  Ap  >  cu2 p3 


(2.27) 


(2.28) 


1 


=>  cos 


(^co  +  £2^  — < 


COs(ftTt) 
-cos  (cut} 


if  z  =  0  Ap  -  cu2 p3 
if  z  =  0  Ap  <  cu2 p3 
if  z  =  0  /\p  >  cu2 p3 


(2.29) 


From  Eq.  (2.29),  and  remembering  the  expression  of  the  eccentricity  provided  in  Eq.  (2.14),  the  expression  of  the 
in-plane  element  /  for  a  zero-displacement  orbit  is 
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f\z=0  -  e  COS  (  +  Q)  =  m  ^ — —  COS  (&7/) 


Considering  now  a  displaced  orbit,  Eq.  (2.11)  gives  the  expression  of  the  RAAN  as 

n|  ={  Q  ,  lfA 

k*°  \2tz  -  Q  if  ~A 

in  which  the  condition 

A  =  (cos  ( mt)  <  0  a  z  >  0)  v  (cos  ( zut )  >  0  a  z  <  0) 

The  expression  of  the  argument  of  perigee,  given  by  Eq.  (2.18),  is 

r  nf  2  if  B 
co  =  { 

[37r/2  if  C 

in  which  the  two  conditions  above  are 

B  =  ^z>0a/j<  m2p2^p2  +z2  j  <  0  a  p>  m2p2^p2  +  z2  j 

< 

C  =  (z>0ap>  m2p2^p2  +z2  jvfz  <  Oa  p<  m2  p2*\p2  +  z  j 

From  Eqs.  (2.31)-(2.33),  the  cosine  of  the  longitude  of  perigee  is 

cos(;r/2  +  Q*)  =  -sin(f2*)  if  A  a  B 
cos(3;r/2  +  Q*)  =  sin(Q*)  ifAAC 

=> 

cos  (5^/2  -  Q*  )  =  sin  (Q*)  H-AaB 
cos  (171/2  -  Q* )  =  -sin(Q*)  if  ~A  aC 

f  esin(Q*)  if  (AaC)v(~Aa5) 
=^>  f\  =  ecos(&>  +  Q)  =  J 

'"u  |-£sin(Q*)  if  (Aa5)v(-AaC) 

After  some  mathematical  manipulations,  Eq.  (2.36)  becomes 

e  cos  ( zut )  sgn  ( z )  if  B 


cos(a>  +  Q)  = 


f\zf0=ecos(co  +  Q)  = 


?cos(^7f)sgn(z)  if  C 


(2.30) 


(2.31) 


(2.32) 


(2.33) 


(2.34) 


(2.35) 


(2.36) 


(2.37) 


Recalling  the  expression  of  the  eccentricity  given  in  Eq.  (2.14)  and  the  conditions  shown  in  Eq.  (2.34),  it  is  possible 
to  demonstrate  that  the  expression  of  the  in-plane  element  /  for  out-of-plane  orbits  shown  in  Eq.  (2.37)  becomes 


,  m2p1*p>rTz2  -n  (  A 

f  L  = — - cosM 

p 


(2.38) 


Note  that  the  expression  for  the  in-plane  element  /  in  case  of  zero -displacement  orbits  shown  in  Eq.  (2.30)  is  a 
particular  case  of  Eq.  (2.38).  Therefore,  the  general  expression  for  the  in-plane  element  /  is  given  by  Eq.  (2.39). 


uj2p2s[prW  -n  (  t 

f  = - - - cos  ( zut ) 


(2.39) 
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3.  In-plane  element  (  g  ) 

As  for  the  case  of  the  in -plane  element  /  ,  let  us  divide  the  derivation  of  the  general  mapping  for  the  in -plane 
element  g  considering  the  in-plane  and  out-of-plane  orbits  separately  at  first.  Starting  by  considering  a  zero-displace¬ 
ment  orbit  and  remembering  the  value  of  the  longitude  of  perigee  shown  in  Eq.  (2.28),  the  value  of  the  sine  of  the 
longitude  of  perigee  is  given  by 


sin  (<z>  +  Q)  = 


0 


if  z  =  0  a  ju  =  m2  p3 


sin (W)  if  z  =  0  a//  <  m2  p2 
-sin (tut)  if  z  =  0  a  ju  >  m2  p2 


(2.40) 


From  Eq.  (2.40),  and  remembering  the  expression  of  the  eccentricity  shown  in  Eq.  (2.14),  the  expression  of  the  in¬ 
plane  element  g  for  a  zero-displacement  orbit  is 


2  3  _ 

g\  _Q  =  esin (&>  +  Q)  =  — ^ — —  sinfW) 


(2.41) 


Following  the  same  arguments  as  those  for  the  derivation  of  Eq.  (2.35),  the  sine  of  the  longitude  of  perigee  is 

sin(;r/2  +  Q*)  =  cos^Q*)  if  A  a  B 
sin  ^3tt/  2  +  Q*  )  =  -cos  (Q*)  ifAAC 


sin  (&>  +  Q)  = 


in(5;r/2-Q*)  =  cos^Q*)  if~A/\B 
sin(7;r/2-Q*)  =  -cos^Q*)  if~AAC 


(2.42) 


sL=esin(®+QH 


cos (Q  j  if  (AaS)v(-AaB) 
(Q*)  if  (AaC)v(~a1aC) 


(2.43) 


-e  cos 


Note  that  both  the  conditions  in  Eq.  (2.43)  do  not  depend  on  the  condition  A  ,  so  that  Eq.  (2.43)  can  be  rewritten  as 


sL=esin(®+Q): 


ecos  (Q*)  if  B 
|^-ecos(Q*)  if  C 


(2.44) 


Recalling  the  definition  of  the  angle  Q*  =  arccos(sin(W)sgn(z)) ,  Eq.  (2.44)  can  be  further  simplified  and  re¬ 
written  as 


=  esm 


(&>+Q)  = 


esin(c7t)sgn(z)  if  B 
-esin(ft7t)sgn(z)  if  C 


(2.45) 


As  for  the  in-plane  element  /  ,  recalling  the  expression  of  the  eccentricity  given  in  Eq.  (2.14)  and  the  conditions 
shown  in  Eq.  (2.34),  it  is  possible  to  demonstrate  that  the  expression  for  the  in-plane  element  g  for  out-of-plane  orbits 
shown  in  Eq.  (2.45)  becomes 


tU +Z~  ~M„:„ 


-sin 


A 


M 


(2.46) 
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Note  that  the  expression  for  the  in-plane  element  g  in  case  of  the  zero-displacement  orbit  shown  in  Eq.  (2.41)  is 
a  particular  case  of  Eq.  (2.46).  Therefore,  the  general  expression  for  the  in-plane  element  g  is  given  by  Eq.  (2.47). 


m^p2^p^  -  p 

P 


sin 


(mt) 


(2.47) 


4.  Out-of-plane  element  (h  ) 

Before  deriving  the  expression  for  the  out-of-plane  element  h  ,  let  us  recall  the  following  identity. 


tan 


fx)_  l-cos(x) 

V  2  J  sin(x) 


(2.48) 


Recalling  the  formulation  of  the  inclination  shown  in  Eq.  (2.8)  as  i  =  arccos 
tan(//2)  is,  therefore,  given  by 


P 


p2  +z2 


,  the  expression  of 


1- 


P 


“".2.  ' 


l 


P2+Z 2 


sin 


f  r 

P 


arccos 


i 


p2  +  z2 


1- 


P 


4 


p2  +z2 


p2  +z2 


(2.49) 


JJ 


tan 


a- 


p~4 

2  +  z2 

\z\ 

(2.50) 


From  Eq.  (2.50),  and  recalling  the  value  of  cos(Q)  =  sin(srf)sgn(z)  introduced  in  Eq.  (2.1 1),  the  expression  of 
the  out-of-plane  element  h  is 


h  =  — 


sm(mt)f-4prZzr'^ 


(2.51) 


Note  that  there  are  no  singularities  in  Eq.  (2.51)  in  the  case  of  zero-displacement  orbits.  It  can  be  seen,  in  fact,  that 

1  i  m  h  =  0  (2.52) 

z— >0 


5.  Out-of-plane  element  (  k  ) 

For  the  expression  of  the  out-of-plane  element  k  ,  the  formulation  of  sin(Q)  is  needed.  From  Eq.  (2.11),  it  is 


sin(Q) 


sin  (Q* )  =  ^1-sin2  (zz?t)sgn2  (z) 
sin^2^--Q*)  =  -^1-sin2  (^)sgn2  (z) 


if  z  =  0  v  (cos (mt )  <  0  a  z  >  0)  v  (cos (mt)  >  0  a  z  <  0) 
otherwise 


(2.53) 
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Defining  the  auxiliary  variable  w  :=  ^1-sin2  (mt)  sgn2  (z)  ,  the  formulation  of  the  out-of-plane  element  k  is 

-w  if  z  =  0  v  (cos  [mt )  <  0  a  z  >  0)  v  (cos  ( mt )  >  0  a  z  <  0) 


k  =  tan(  ^  Jsin(Q)  =  < 


p-4 

rj- 

2  +  z2 

z 

p-4 

2  ,  2 

+  Z 

|z| 

(2.54) 


w  otherwise 


Let  us  study  the  two  branches  of  the  solution  separately  starting  from  the  first  one.  If  the  first  condition  in  Eq. 
(2.54)  is  verified,  there  is 


p-4 

n 

2  ,  2 

+  Z 

|z 

^1-sin2  (ftTt)sgn2  (z)  = 


ifz  =  0 


\cos(m)\(P-47T7) 


if  z  >  0 


|cos  {jut )|  | p  -  *Jp2  +z2  j 


if  z  <  0 


(2.55) 


k  - 


cos 


0  if  z  =  0 

{&t)[p-Jp2  +  z2} 


otherwise 


On  the  contrary,  if  the  second  condition  of  Eq.  (2.54)  is  verified  then 

i(m)\(p-'[p2 +£2) 


k  = 


p-4 

r 

2  +  z2 

z 

if  z  >  0 


cos 


if  z  <  0 


cos 


k  =  - 


(mt)^p-^p1  +z2^ 


Therefore,  Eq.  (2.54)  can  be  rewritten,  by  unifying  Eq.  (2.56)  and  Eq.  (2.58),  as 

0  if  z  =  0 

cos(mt)^p-^[p2  +z2  j 


k  = 


otherwise 


(2.56) 


(2.57) 


(2.58) 


(2.59) 
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However,  it  can  be  demonstrated  that 


lim 

z— >0 


COS 


(znt)[^p - p2  +z2 ) 


=  0 


(2.60) 


Therefore,  the  first  condition  in  Eq.  (2.59)  is  actually  redundant  and  unnecessary.  The  general  expression  for  the 
out-of-plane  element  k  is  therefore  given  by 


cos 


k  =  - 


(2.61) 


6.  True  longitude  ( L  ) 

In  order  to  derive  the  general  formulation  for  the  true  longitude  L  ,  let  us  consider  three  cases  z  =  0 ,  z  >  0 ,  and 
z  <  0 ,  respectively. 

In  the  case  of  zero-displacement  orbits  (i.e.  z  =  0 ),  the  expression  of  the  true  longitude  is  simply  given  by  the 
formulations  of  RAAN,  argument  of  perigee,  and  true  anomaly  given  in  Eq.  (2.23),  as 

l\z0  =m  (2.62) 

In  the  case  of  positive  out-of-plane  displacement  (i.e.  z  >  0 ),  the  formulations  of  the  angles  Cl  ,  co  ,  and  3  are 
the  following. 


■  |Q*  ifcos(W)<0 

z>0  [2 n  -  Cl*  otherwise 

,  U/2  if  p<zu2p2Jp2+z2 

co\  n  =  \  _ _ 

[3 7r/ 2  if  p  >  TU2 p2  y] p 2  +Z2 

=  |o  ii^i<m2p2^prW 

\x  if  p>  m2 p2  sf/TTz2 

From  Eqs.  (2.64)-(2.65)  and  considering  the  angles  lying  in  the  interval  [0,  In)  then 


(2.63) 


(2.64) 


(2.65) 


[®  +  '9Lo=f 

Therefore,  the  expression  of  the  true  longitude  for  positive  out-of-plane  displacement  is 


(2.66) 


L\  =  — +  Q|  = 

lz>0  ^ 


— +  Q*  ifcos(s7r)<0 
2  v  ; 


—  n  -Cl  otherwise 
2 


(2.67) 


In  the  case  of  negative  out-of-plane  displacement  (i.e.  z  <  0 ),  the  formulations  of  the  angles  Cl ,  co  ,  and  3  are 
the  following. 
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°L  - 


co\  n  =  < 

h<0 


Jq* 

if  cos(W)  >  0 

{2^- 

-  Q*  otherwise 

7r/  2 

if  p  >  m2p2^p2  +z2 

3tt/  2 

if  //  <  GT2  p1  p2  +  Z2 

4 


[o  if  p  <  m2p2^jp2  +z2 
[tt  if  /j>m2p2^p2  +z2 


(2.68) 


(2.69) 


(2.70) 


From  Eqs.  (2.69)-(2.71)  and  considering  the  angles  lying  in  the  interval  [0,2;r)  then 

[CD  +  S\<o=\7r 

Therefore,  the  expression  of  the  true  longitude  for  negative  out-of-plane  displacement  is 


(2.71) 


—  71  + Cl'  ifcOs(W)>0 
2  V  7 


-7T-Q  otherwise 
2 


(2.72) 


The  general  expression  for  the  true  longitude  L  is  given  by  unifying  Eq.  (2.62),  Eq.  (2.67),  and  Eq.  (2.72)  as 

mt  if  z  =  0 

if  z  >  0  acos(gj t)  <  0 


L  - 


-  +  Q* 
2 


(2.73) 


-7T-Q *  if  z  >  0  acos (mt)  >  0 
2  v  7 

3  * 

—  71  +  0*  if  z  <  0  ACOs(ftTt)  >  0 
2  v  7 

7  * 

-7T-Q*  if  z  <  0  ACOs(ftTt)  <  0 


A  further  simplification  of  the  general  formulation  of  the  true  longitude  shown  in  Eq.  (2.73)  can  be  carried  out  by 
studying  the  variation  of  the  angle  Q  .  In  order  to  further  simplify  the  expression  of  the  true  longitude,  let  us  start 
recalling  the  formulation  of  the  angle  Q*  =  arccos(sin(W)sgn(z)) .  A  study  on  the  characteristics  of  this  function 

can  be  done  by  plotting  the  value  of  Q  against  the  angle  x  =  mt .  From  Eq.  (2.73)  it  is  evident  that  there  is  no  interest 
in  studying  the  value  Q  .  Therefore,  two  separate  plots  are  shown  in  Figure  2,  considering  the  case  of  positive  and 

lz=0 

negative  vertical  displacement,  respectively. 
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Figure  2.  Plot  of  arccos(sin(x)sgn(z)). 

From  the  plots  shown  in  Figure  2,  it  can  be  verified  that  the  expression  of  Q  in  Eq.  (2.74)  is  true. 


-n/l+m 
5^/2  -mt 
^nfl-mt 
— 3tt/  2  +  mt 


if  z  >  OAcos(mt)  <  0 
if  z  >  0  a  cos  ( mt )  >  0 
if  z  <  0  a  cos  ( mt )  <  0 
if  z  <  0  a  cos  (mt )  >  0 


Therefore,  Eq.  (2.73)  can  be  rewritten  in  a  simpler  form  as 

L  =  mt 


(2.74) 


(2.75) 


Summary 

The  forward  map  from  NKOs  to  osculating  modified  equinoctial  elements  described  in  Section  2. 1.1. 2  is  summa¬ 
rized  as: 


P  = 

f  = 

8  = 


P 

m2p2yjp2  +z2  -ju 
P 

&2P2^P2  +z2  -p 

P 


cos  (mt) 
sin  (mt) 


sm(mt)(p-p>  +  z 
h  = - - - 

z 

cos  (mt ) (p  - 7^+T 
k  = - i 

z 

L  =  mt 


(2.76) 
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2.I.I.3.  Augmented  Integrals  of  Motion 

The  augmented  integrals  of  motion  are  defined  as 

h  Orbital  angular  momentum  vector 

<  e  Eccentricity  vector  (2.77) 

L  True  longitude 


Note  that  there  exists  a  relationship  between  the  two  integrals  of  motion  (i.e.  h  e  =  0 ),  which  reduces  the  number  of 
independent  elements  for  the  orbit  description  from  7  to  6. 

Starting  from  the  expressions  of  the  Cartesian  position  and  velocity  vectors  (Eqs.  (2.5)-(2.6)),  the  mapping  from 
NKO  elements  [z,p,ftjf  to  osculating  augmented  integrals  of  motion  is  shown  in  Eq.  (2.78).  The  formulations  of 

angular  momentum  vector,  eccentricity  vector,  and  true  longitude  are  given  from  Eqs.  (2.7),  (2.15),  and  (2.75),  re¬ 
spectively. 


h  - 


e  - 


-TUp Z  COS  (ftTt) 

-mpz  sin  {tut) 

2 


mp 


_2  2 

zu  p 


p2  +z2 


pCOs(ftTt) 

psin(mt) 

z 


L  =  mt 


(2.78) 


2.1.2.  Sensitivity  Analysis 

A  sensitivity  analysis  has  been  performed  in  order  to  study  the  behavior  of  the  mappings  in  the  presence  of  errors 
in  the  NKO  elements  [z,p,zjf  .  Both  analytical  and  numerical  sensitivity  analyses  have  been  carried  out  to  confirm 
the  results  found. 

The  analytical  sensitivity  analysis  is  centered  on  the  computation  of  the  Jacobian  J  of  the  mapping.  Defining  the 

non-Keplerian  orbit  to  osculating  orbital  elements  map  as  O  and  the  set  of  NKO  elements  y  =  [z,p,tfff  ,  the  analyt¬ 
ical  sensitivity  analysis  is  performed  by  means  of  the  Jacobian  of  the  system  such  that 


dO 

dy 


(2.79) 


In  the  following  sub-sections,  the  definition  of  the  parameters  used  for  the  sensitivity  analysis  and  a  definition  of 
the  structure  of  the  Jacobian  are  given  for  each  set  of  osculating  orbital  elements. 
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2.I.2.I.  Classical  Keplerian  Elements 

The  mapping  ®  from  non-Keplerian  orbit  to  osculating  classical  Keplerian  elements  is  shown  in  Eq.  (2.23).  The 
Jacobian  of  this  mapping  is 


J  KEP 


da 

de 

di 

ao 

dco 

89  ~ 

dz 

dz 

dz 

dz 

dz 

dz 

da 

de 

di 

5Q 

dco 

89 

dp 

dp 

dp 

dp 

dp 

dp 

da 

de 

di 

50 

dco 

89 

dzu 

dzu 

dzu 

dzu 

dzu 

dzu 

(2.80) 


The  elements  of  the  Jacobian  are  not  shown  here  for  the  sake  of  conciseness. 

It  is  important  to  note  that  some  of  the  formulations  found  for  the  NKO  to  osculating  classical  Keplerian  elements 
map  are  expressed  through  a  piecewise  law,  which  creates  difficulties  in  dealing  with  analytical  sensitivity  analysis. 
It  is  not  possible,  for  instance,  to  analytically  differentiate  the  expression  of  right  ascension  of  the  ascending  node, 
argument  of  perigee,  and  true  anomaly.  Therefore,  for  what  concerns  these  elements,  only  a  study  based  on  numerical 
sensitivity  analysis  has  been  carried  out. 

The  parameters  used  to  test  the  sensitivity  of  the  mapping  to  uncertainties  in  the  NKO  parameter  y.  are  as  follows. 

•  Relative  error  of  the  semi -major  axis  due  to  uncertainties  in  the  NKO  elements. 

\a-a\  \a-(a  +  Sa)\  \Sa\  \ja,yt  '8y\ 


8a,yt  =- 


(2.81) 


Absolute  error  of  the  eccentricity  due  to  uncertainties  in  the  NKO  elements. 

£e,y.  :=\e-e\=\e-(e  +  8e)\  =  \8e\=\jeyi-8yi\  (2.82) 

Normalized  error  of  the  angles  £  =  j i,  O,  co,  9\  due  to  uncertainties  in  the  NKO  elements. 


8$. 


14*  •‘ty 


£Z,yt  ~  ' 

n  n 


(2.83) 


2.I.2.2.  Modified  Equinoctial  Elements 

The  mapping  ®  from  non-Keplerian  orbit  to  osculating  modified  equinoctial  elements  is  shown  in  Eq.  (2.76).  The 
Jacobian  of  this  mapping  is 


~8p 

df 

dg 

8h 

dk 

dL  " 

dz 

dz 

dz 

dz 

dz 

dz 

dp 

df 

dg 

8h 

dk 

dL 

dp 

dp 

dp 

dp 

dp 

dp 

dp 

df 

dg 

dh 

dk 

dL 

dzu 

dzu 

dzu 

dzu 

dzu 

dzu 

(2.84) 


Again,  the  elements  of  the  Jacobian  are  not  shown  here  for  the  sake  of  conciseness. 
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The  parameters  used  to  test  the  sensitivity  of  the  mapping  to  uncertainties  in  the  NKO  parameter  yt  are  as  follows. 
•  Relative  error  of  the  semi-latus  rectum  due  to  uncertainties  in  the  NKO  elements. 


\p~p\  _  \P~{P  +  SP)\  _  \Sp\  _  \Jp,y,  -^1 

P  P  ~~~  P 

Absolute  error  of  the  elements  x  =  {/,  g,h,k }  due  to  uncertainties  in  the  NKO  elements. 

£X  Ji  :=  |x-x|  =  |x-(x  +  £x)|  =  \Sx\  =  |  Jx  y_  -Syi  | 

Normalized  error  of  the  true  longitude  L  due  to  uncertainties  in  the  NKO  elements. 

sl  K,  M. 


£L,y,  ■= 


(2.85) 


(2.86) 


n 


n 


(2.87) 


2.I.2.3.  Augmented  Integrals  of  Motion 

The  mapping  <D  from  non-Keplerian  orbit  to  osculating  augmented  integrals  of  motion  is  shown  in  Eq.  (2.78). 
The  Jacobian  of  this  mapping  is 


dhx 

8hy 

dhz 

dex 

dey 

dez 

dL  ~ 

dz 

dz 

dz 

dz 

dz 

dz 

dz 

dhx 

dhy 

dhz 

dex 

dey 

dez 

dL 

dp 

dp 

dp 

dp 

dp 

dp 

dp 

dhx 

dhy 

dhz 

dex 

dey 

dez 

dL 

dm 

dm 

dm 

dm 

dm 

dm 

dm 

(2.88) 


Again,  the  elements  of  the  Jacobian  are  not  shown  here  for  the  sake  of  conciseness. 

The  parameters  used  to  test  the  sensitivity  of  the  mapping  to  uncertainties  in  the  NKO  parameter  y.  are  as  follows. 
•  Magnitude  of  the  difference  vector  between  the  nominal  value  of  the  angular  momentum  and  that  related 
to  uncertainties  in  the  NKO  elements.  Relative  error. 


||*-*||  ||ft -(A+£A)||  \\sh\\  \jh,yi  'Syt | 


(2.89) 


h'y‘  h  h  h  h 

Magnitude  of  the  difference  vector  between  the  nominal  value  of  the  eccentricity  vector  and  that  related 

to  uncertainties  in  the  NKO  elements.  Absolute  error. 

£e,yi  '■=  \\e-£\\  =  Ik -(e  +  ^)ll=  IW  HK*  (2-90) 

Normalized  error  of  the  true  longitude  L  due  to  uncertainties  in  the  NKO  elements. 


^  := 


SL 


J 


L,yt 


n 


n 


(2.91) 
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2.1.3.  Numerical  Test  Cases 

Two  different  sets  of  test  cases  have  been  considered.  The  first  set  has  been  chosen  to  validate  the  forward  map¬ 
pings  and  to  study  the  behavior  of  the  osculating  elements.  The  parameters  have  been  selected  for  ease  of  illustration. 
Therefore,  the  acceleration  required  to  generate  the  NKOs  considered  is  not  constrained  and  so  can  be  well  above  the 
technical  limitations  of  the  current  low-thrust  technology.  The  second  set  has  been  chosen  to  study  the  impact  of 
uncertainties  in  the  non-Keplerian  orbit  properties  on  the  classical  orbital  elements.  These  two  sets  are  described  in 
Section  2. 1.3.1  and  Section  2. 1.3. 2,  respectively. 


2.I.3.I.  Mappings 

15  different  scenarios  have  been  tested  in  order  to  both  verify  the  validity  of  the  mappings  derived  in  the  previous 
sections  and  to  study  the  behavior  of  the  osculating  elements,  as  detailed  below.  All  results  of  the  test  cases  have  been 
numerically  verified  by  means  of  MATLAB. 

In  the  following  sub-section,  all  test  cases  are  described.  Figure  3  shows  the  orbits  of  all  these  test  cases  used  in 
this  study. 


Figure  3.  Earth-centered  view  of  the  orbits  of  all  test  cases  used. 
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Test  cases 


1.  GEO 

The  first  test  case  is  a  conventional  Geostationary  Earth  Orbit  (GEO)  described  by 


GEO: 


^geo  0 

Pgeo  —  rGEo  ~  42, 1 64  km 


(2.92) 


1  rGEO 


2.  Type  1  NKO 

Two  test  cases  are  considered  for  the  Type  1  NKO. 

Type  1  NKO : 


h  =  rGEo  sin(M) 
P\  =  rGE0  cos(A2) 
1  =  '\p/VGEO 


(2.93) 


67, 


The  term  A/L  is  the  angle  between  the  position  vector  r  and  the  X —Y  plane.  In  this  case,  two  values  are  con¬ 
sidered  so  that  both  positive  and  negative  out-of-plane  displaced  NKOs  are  taken  into  account. 


AA  =  ±20  deg 


z1  =  ±0.34rGEO  =  ±14,000  km 
px  -  0.94 rGEO  =  39,600  km 


(2.94) 


3.  Type  2  NKO 

Four  test  cases  are  considered  for  the  Type  2  NKO. 


Type  2  NKO : 


z2  =  {±5-/?0,±lO*/?0} 


Pi  =  '* 


(2.95) 


67 


=TTr 


4.  Type  3  NKO 

Eight  test  cases  are  considered  for  the  Type  3  NKO.  Six  of  them  are  characterized  by  an  out-of-plane  displacement 
while  the  other  two  are  in-plane  displaced  NKOs. 


Type  3  NKO 

(out-of-plane  displacement) 


Type  3  NKO 
(in-plane  displacement) 


^3 op  ~  ~rGEO 

<  P3op  =  {0-9,  1?  l'l}‘rG£0 
m3op  =  '\p/rGEO 

hip  =0 

P3ip  ~  {^-9,  1.1  }‘rGEO 

m3ip  ~  ' \p/rGEO 


(2.96) 


(2.97) 
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Results  and  discussion 

Figure  4  shows  the  Type  1  NKO  with  >°  described  by  Eq.(2.93),  together  with  its  osculating  Keplerian  orbits, 

corresponding  to  the  instantaneous  orbital  elements  of  the  spacecraft  if  the  thrust  is  nulled.  In  this  case,  it  is  shown 
how  the  semi-major  axis  and  eccentricity  of  the  osculating  orbits  do  not  change,  but  the  osculating  orbit  plane  rotates 
around  the  vertical  axis.  The  spacecraft,  being  always  at  the  apogee  of  the  osculating  Keplerian  orbits,  describes  the 
desired  NKO.  Therefore,  the  envelope  resulting  from  the  osculating  Keplerian  orbits  after  an  entire  orbit  is  well  ap¬ 
proximated  by  a  truncated  cone,  characterized  by  height  H  and  radii  of  the  bases  rv  r2 ,  as  follows. 

rx  -  a(\-e)cosi 

<r2=p  (2.98) 

H  -  |z|  +  |a(l-e)sin/| 


Figure  4.  Type  1  out-of-plane  displaced  GEO  (solid  black  line)  and  osculating  Keplerian  orbits.  In  gray  is  shown  the 
truncated  cone  envelope  resulting  from  the  osculating  Keplerian  orbits  after  an  entire  orbit. 


In  the  following  sub-sections,  the  results  of  the  mappings  using  each  set  of  osculating  elements  are  shown  and 
discussed. 
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Table  1  shows  the  summary  of  the  results  of  the  non-Keplerian  orbit  to  classical  Keplerian  elements  map  for  all 
the  test  cases  taken  into  account. 


Table  1.  Summary  of  the  results  of  the  non-Keplerian  orbit  to  classical  Keplerian  elements  map.* 


GEO 

Type  1  NKO 

Type  2  NKO 

Type  3  NKO 
(out  of  plane) 

Type  3  NKO 
(in  plane) 

a  [  rGEO  ] 

1 

0.9 

~10  (z2=±10/?e) 

~1.7  (z2=±57?e) 

~1.5  (Aop 
-  2.4  ( p3op  —  rG£o) 
-7.4  [P3op  —  l-lrG£6>) 

°-7  {p3ip=0.9rGEO) 

1 .6  (P}ip  =  1  •  1  rGEO  ) 

e 

0 

0.12 

0.8  (z2=±10fle) 

0.25  (z2  =  ±5/?e) 

0.09  i^P^op  —  ®'9rGEO  ) 
0.4  i^P^op  ~  rGEO  ) 

0*8  (P3op  =1-Keo) 

0-27  (p3ip  =  0.9igfo  ) 
0.33  (/7?ip  =  1-1?^£0) 

i  [deg] 

0 

20 

56  (z2  =  ±107?e) 

37  (z2  =  ±5/?e) 

48  ( P3op  ~  ®’9rGEO  ) 

45  ( P3ap  =  rGEO  ) 

42  (  p3op  ~  1  1 rGEO  ) 

0 

O  [deg] 

0 

90  Z  (z,  <  0) 
270  Z  (z,  >  0) 

90  Z  (z2  <  0) 

270  Z  (z2  >  0) 

90  Z  ( ziop  <  0) 

270  Z  (z3op  >  0) 

0 

co  [deg] 

0 

90  (zj  <  0) 
270  (-,  >0) 

270  (z2  <  0) 

90  (z2  >  0) 

270  (z3op<0) 

90  (z3op  >  0) 

180  Z  ( P3ip<rGE0 ) 

0  Z  ( Pyp  >  rGEO  ) 

9  [deg] 

OZ 

180 

0 

0  (P3op  ~  rGEO  ) 

180  (p3op  <  rGEO  ) 

180  (/?3*p  <  rGEO  ) 

^  {P3ip  >  rGEO  ) 

*  The  symbol  is  used  with  the  meaning  “linearly  increasing  with  slope  m  The  entry  T  ,  for  instance,  means  that  the  value  of  the  generic 
angle  z  in  time  is  T  (f )  =  To  +  mt  . 
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2.  Modified  equinoctial  elements 

Table  2  shows  the  summary  of  the  results  of  the  non-Keplerian  orbit  to  modified  equinoctial  elements  map  for  all 
the  test  cases  taken  into  account. 


Table  2.  Summary  of  the  results  of  the  non-Keplerian  orbit  to  modified  equinoctial  elements  map.* 


GEO 

Type  1  NKO 

Type  2  NKO 

Type  3  NKO 
(out  of  plane) 

Type  3  NKO 
(in  plane) 

p 

\_rGEO  ] 

1 

0.88 

~  3.3  (z2=±10 R9) 
~1.6  (z2  =±5 Re) 

~  1-5  {P3op  =  ®'9rGEO  ) 

~  ^  ( p3op  =  rGEO  ) 

~  2.7  (p3op  =1-1  rGEo) 

(P3ip  ~  ®'9rGEO  ) 
l-^  {P3ip  =1*1  rGEo) 

/ 

0 

f  =  f0  cos  (m) 
with  f0<  0 

f=f0cos(m)  with  f0>  0 

f  =  f0  cos  (W)  with 

\fo  <  ^  {Plip  <  rGEO  ) 

\fo  >  ^  {Plip  >  VGEO  ) 

8 

0 

g  =  So  sin(W) 

with 

f  P  ^ 

g  t<^~  <0 

9 

V  z  2 

/  rji  \ 

g  =g0  sin  (mt)  with  g  t<  >0 

V  Z  ) 

8  =  g0sin(m)  with 

< 

^<Zf)>0  ( P3ip>rGEO ) 

h 

0 

/z  = /Zq  sin(W)  with  < 

h{t<~ ( z< °) 

h[t<~ y]>0  (z>0) 

0 

k 

0 

(fc0>0  (z<0) 

k  =  kn  cost  tut)  with  {  ;  ; 

0  V  ;  \k0<0  (z>0) 

0 

L 

[deg] 

o/1 

0^ 

O/1 

O/1 

o/1 

*  The  symbol  is  used  with  the  meaning  “linearly  increasing  with  slope  w  The  entry  T(]  P  ,  for  instance,  means  that  the  value  of  the  generic 
angle  z  in  time  is  T  (t )  =  To  +  mt  . 


Figures  5-6  show  the  evolution  over  one  orbit  of  the  in -plane  modified  equinoctial  element  /  and  the  out-of- 

plane  modified  equinoctial  element  h ,  respectively.  It  can  be  seen  that  the  amplitudes  of  the  oscillations  depend 
respectively  on  the  eccentricity  and  inclination  of  the  osculating  orbits,  as  defined  in  Eq.  (2.24). 
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Figure  5.  Evolution  of  the  osculating  in-plane  modi-  Figure  6.  Evolution  of  the  osculating  out-of-plane 

fied  equinoctial  element/ over  one  orbit.  modified  equinoctial  element  h  over  one  orbit. 


3.  Augmented  integrals  of  motion 

Figure  7  shows  the  cone  created  by  the  osculating  angular  momenta  due  to  the  rotation  of  the  osculating  Keplerian 
orbits.  The  height  of  the  cone  and  its  base  radius  depend  on  both  vertical  displacement  z  and  radius  of  the  orbit  p  . 
For  example,  the  two  Type  2  NKOs  with  the  largest  vertical  displacement  are  those  with  the  largest  base  radius.  On 
the  other  hand,  the  two  Type  2  NKOs  with  \z2  |  =  5  •  R@  have  a  smaller  base  radius  than  the  out-of-plane  Type  3  NKOs, 

which  have  | z3op  |  =  rGEO .  Among  the  Type  3  NKOs,  the  orbits  with  the  lowest  radius  p3op  have  the  largest  base  radius, 

followed  by  the  orbits  with  the  largest  radius  and,  last,  those  orbits  with  p3op  =  z3op  .  The  osculating  angular  momenta 

for  the  in-plane  Type  3  NKOs  is  parallel  to  the  angular  momentum  of  GEO  and  therefore  not  clearly  distinguishable 
in  the  plot. 


GEO 

-  -‘■Type  1  NKO  -  Displaced  GEO 

Type  2  NKO  -  Synchronous  with  a  GEO 
Type  3  NKO  -  Displaced  geosynchronous  (out  of  plane) 
Type  3  NKO  -  Displaced  geosynchronous  (in  plane) 


.2  Y  [km  /s] 


X  [knrr/s] 


Figure  7.  Time  evolution  of  the  osculating  angular  momentum  vector.  All  the  test  cases  studied  are  shown. 
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Figure  8  shows  the  double  cone  created  by  the  osculating  eccentricity  vectors  due  to  the  rotation  of  the  osculating 
Keplerian  orbits.  The  height  of  the  cone  and  its  base  radius  depend  on  both  the  vertical  displacement  z  and  the  radius 
of  the  orbit  p  .  The  osculating  eccentricity  vector  depend  on  the  sign  of  the  vertical  displacement,  whereas  the  oscu¬ 
lating  angular  momentum  does  not.  Moreover,  in  this  case,  the  behavior  is  different  with  respect  to  the  angular  mo¬ 
mentum,  as  can  be  better  seen  from  Figure  9.  For  example,  the  two  Type  2  NKOs  with  the  largest  vertical  displacement 
are  not  those  with  the  largest  base  radius,  but  are  those  with  the  largest  value  of  eccentricity.  On  the  other  hand,  the 
two  Type  3  NKOs  with  p3op  =1.1  -rGE0  have  the  largest  base  radius.  The  osculating  eccentricity  vectors  for  the  in- 


plane  Type  3  NKOs  lie  on  the  X—Y  plane,  while  the 

eGEO=[l  0  Of- 


Figure  8.  Time  evolution  of  the  osculating  eccentricity 
vector.  All  test  cases  studied  are  shown. 


eccentricity  vector  in  the  GEO  case  is  arbitrarily  set  to 


Figure  9.  Time  evolution  of  the  osculating  eccentricity 
vector.  X-Y  view. 


2.I.3.2.  Sensitivity  Analysis 

Three  test  cases  have  been  used  to  investigate  the  sensitivity  of  the  mappings  to  a  change  in  the  NKO  elements 
[z,  p,  z&f  .  All  the  test  cases  are  chosen  as  feasible  mission  scenarios,  with  reasonable  required  low-thrust  acceleration, 
and  are  all  displaced  GEO. 

The  first  test  case  is  an  out-of-plane  displaced  GEO  Type  1  NKO,  defined  as  in  Eq.  (2.93).  The  displacement  angle 
AA  has  been  chosen  such  that  the  spacecraft  hovers  above  GEO  at  a  vertical  displacement  equal  to  twice  the  station - 
keeping  box.  Station-keeping  regulations  currently  require  a  station -keeping  box  height  of  A  A  e  [0.05,  O.ljdeg  .4  The 
worst-case  scenario  is  taken  into  account  here  so  that  Eq.  (2.93)  becomes 


Out-of-plane  displaced  GEO 
(Type  1  NKO) 


z1  =  rGEO  sin(0.2deg)  =  147  km 
'  P\  —  rGEo  cos(0.2deg)  ~  rGEO 
®i 


(2.99) 


A.  Peloni^|^T^^p,(^^i^proved  for  pub|jc  re|ease.  distribution  unlimited. 


Page  26 


Osculating  Keplerian  Elements  for  Highly  Non-Keplerian  Orbits  -  Final  technical  report 


The  second  test  case  is  an  in-plane  displaced  GEO  Type  3  NKO.  A  positive  in-plane  displacement  of 
Ap3  =  147  km  is  considered.  This  value  has  been  chosen  to  be  the  same  as  the  out-of-plane  displacement,  as  in 
Heiligers  at  al.4 

z3  =  0 

p3  =  rGE0  +147  km  (2.100) 


In-plane  displaced  GEO 
(Type  3  NKO) 


The  last  test  case  has  been  chosen  in  order  to  verify  the  behavior  of  the  map  in  case  of  very  small  vertical  displace¬ 
ments,  such  that  a  small  error  in  the  out-of-plane  displacement  z1(2)  causes  the  orbit  to  cross  the  z  =  0  plane.  Such  a 

test  case  has  been  chosen  as  a  feasible  mission  scenario  for  a  displaced  GEO  Type  1  NKO  with  a  small  vertical 
displacement,  defined  in  Eq.  (2.93).  The  displacement  angle  AA  has  been  chosen  so  that  the  actual  out-of-plane 
displacement  is  zl(2)  ~  101  km  .  Therefore,  Eq.  (2.93)  becomes 


Out-of-plane  displaced  GEO 
(Type  1  NKO) 


Zi(2)  =  rGEO  sin  (-10  4  deg)  =  -0.08  km 
<  A(2)  =  rGEo  cos(— 10  deg)  ~  rGEO 

mi(2)  =  *\p/rGEO 


The  relative  error  of  the  NKO  elements  is  considered  fixed  and  constant,  as  follows. 

JlO  3  (0.1%)  First  two  test  cases 

<  z  [2  (200%)  Third  test  case 

8 p  ~  10  3  (0.1%) 


(2.101) 


(2.102) 


Note  that  the  relative  error  related  to  the  out-of-plane  displacement  z1(2)  for  the  third  test  case  is  much  larger  than 
that  considered  for  the  other  two.  This  has  been  done  in  order  to  have  a  crossing  of  the  z  =  0  plane.  In  fact,  a  200% 
error  on  z1{2)  for  the  third  test  case  corresponds  to  the  same  displacement  as  the  original  but  with  a  different  sign,  as 
shown  in  Eq.  (2.105). 

The  values  of  the  absolute  errors  in  the  three  test  cases  are  as  follows. 


Out-of-plane  displaced  GEO 
(Type  1  NKO) 


8zx  =  £zzx  ~  0.15  km 
<  8px  =  sppx  «  42.2  km 
8m]  -£UJmx  «  0.36  deg/day 


In-plane  displaced  GEO 
(Type  3  NKO) 


Sz3  =  £zzx  «  0.15  km 
<  Sp3  =  s pp3  «  42.3  km 
8m 3  =smtu 3  ~  0.36  deg/day 


Out-of-plane  small-displacement  GEO 
(Type  1  NKO) 


8zl(2)  =  8Z  |z1(2)  |  «  0.16  km 

<  +K2,  =  £pPi( 2)  *  42-2  km 

=  £V&1(2)  ~0-36  deg/day 


(2.103) 


(2.104) 


(2.105) 
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Note  that  the  term  Sz3  in  Eq.  (2.104)  depends  on  the  displacement  zx  of  the  out-of-plane  displaced  GEO.  This  is 
made  on  purpose  in  order  to  have  a  non-zero  error  for  the  out-of-plane  displacement. 

The  analytical  sensitivity  analysis  has  been  described  in  Section  2.1.2.  The  numerical  sensitivity  analysis  has  been 
performed  by  computing  the  values  of  the  osculating  elements  from  the  perturbed  NKO  elements  and  then  computing 
the  errors  as  described  for  the  analytical  case. 


Results  and  discussion 

1.  Classical  Keplerian  elements 

The  sensitivity  analysis  study  demonstrates  that  the  mapping  that  uses  the  classical  Keplerian  elements  is  ex¬ 
tremely  sensitive  to  uncertainties  in  the  NKO  elements.  This  is  mainly  due  to  the  piecewise  formulation  of  the  map 
itself.  A  small  error  on  either  p  or  m  ,  for  example,  can  cause  the  value  of  the  argument  of  perigee  co  to  jump  from 
one  branch  to  another  (see  Eq.  (2.23)).  Moreover,  only  a  numerical  sensitivity  analysis  can  be  carried  out  on  half  of 
the  classical  Keplerian  elements. 

For  the  sake  of  conciseness,  only  the  most  interesting  plots  of  the  sensitivity  analysis  study  are  shown  here.  Table 
3  shows  a  summary  of  the  results  of  the  sensitivity  analysis  study  in  the  case  of  classical  Keplerian  elements. 


Table  3.  Summary  of  the  results  of  the  sensitivity  analysis  study  in  case  of  classical  Keplerian  elements.* 


Errors  due  to 

uncertainties 

Test 

case 

Sz 

Sp 

5vj 

all 

<  10“6% 

0.4% 

0.2% 

all 

<10“8 

3xlO“3 

2xl0“3 

£i 

1 

i  <r4% 

^t 

1 

o 

0 

2 

~o 

3 

0 

0 

1 

0 

0 

linearly  increasing  up  to  0.2% 

2 

up  to  100% 

0 

3 

100% 

linearly  increasing  up  to  0.2% 

1 

0 

100% 

100% 

2 

up  to  100% 

0 

linearly  increasing  up  to  0.2% 

3 

100% 

100% 

100% 

1-3 

0 

100% 

100% 

2 

0 

0 

Note  that  sometimes  ~0  is  used  against  0.  This  is  done  to  stress  the  difference  between  a  null  value  and  a  very  small  (numerical)  value. 


A. 


Peloni^|^Tj^^p|(^^i^pr0ved  for  pubMc  release:  distribution  unlimited. 


Page  28 


Osculating  Keplerian  Elements  for  Highly  Non-Keplerian  Orbits  -  Final  technical  report 


Figure  10  shows  the  evolution  of  the  osculating  RAAN  Q  over  one  orbit.  The  test  case  for  the  small  displacement 
GEO  (test  case  3)  is  shown.  The  nominal  value  and  the  values  due  to  the  uncertainties  on  the  NKO  elements  are 
considered.  It  is  shown  that  an  error  in  the  out-plane  displacement  causes  the  RAAN  to  have  a  constant  180-deg  error 
due  to  the  piecewise  formulation  of  the  map  (Eq.  (2.23)).  Figure  1 1  shows  a  zoom  on  the  Y  axis  of  the  normalized 
error  of  the  RAAN  for  the  same  test  case.  It  is  shown  how  an  error  on  the  orbital  angular  velocity  causes  the  RAAN 
(and,  therefore,  the  orbit  plane)  to  rotate  at  a  different  velocity.  The  first  test  case  shows  the  same  behavior.  Similarly, 
the  argument  of  perigee  shows  the  same  behavior  in  case  of  planar  orbit  (i.e.  test  case  2). 


Figure  10.  Evolution  of  the  osculating  RAAN  over 
one  orbit.  Very  small  displacement  GEO  (test  case  3). 


Figure  11.  Normalized  error  of  the  RAAN  for  the 
case  of  very  small  displacement  GEO  (test  case  3).  Zoom 
on  Y  axis. 
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2.  Modified  equinoctial  elements 

For  what  concerns  the  modified  equinoctial  elements,  the  sensitivity  analysis  study  demonstrates  that  the  mapping 
that  uses  these  elements  is  very  robust  to  uncertainties  in  the  NKO  elements. 

For  the  sake  of  conciseness,  only  the  most  interesting  plots  of  the  sensitivity  analysis  study  are  shown  here.  Table 
4  shows  a  summary  of  the  results  of  the  sensitivity  analysis  study  in  the  case  of  modified  equinoctial  elements. 


Table  4.  Summary  of  the  results  of  the  sensitivity  analysis  study  in  case  of  modified  equinoctial  elements.* 


Errors  due  to 

Test 

Sz 

Sp 

8m 

uncertainties 

case 

all 

~0 

0.4% 

0.2% 

£f 

all 

~0 

<3xl0“3 

<2xl0“3 

Ss 

all 

~0 

<3xicr3 

<2x10  3 

1 

<2xl0“6 

<10“5 

£h 

2 

<2xl0“6 

0 

0 

3 

<10“9 

<6xl0“9 

1 

<2xl0“6 

<9xl0“6 

£k 

2 

<2xl0“6 

0 

0 

3 

<10“9 

<5xl0“9 

£l 

all 

0 

0 

linearly  increasing  up  to  0.2% 

Note  that  sometimes  ~0  is  used  against  0.  This  is  done  to  stress  the  difference  between  a  null  value  and  a  very  small  (numerical)  value. 


Figure  12  shows  the  evolution  over  one  orbit  of  the  in -plane  modified  equinoctial  element  /  .  The  test  case  for 
the  out-of-plane  displaced  GEO  (test  case  1)  is  shown.  The  nominal  value  and  the  values  due  to  errors  in  the  NKO 
elements  are  shown.  It  can  be  seen  that  the  error  on  the  in-plane  element  /  due  to  uncertainties  in  the  NKO  elements 
is  of  the  same  order  of  magnitude  as  the  uncertainties  themselves.  The  same  behavior  is  noted  for  the  other  in-plane 
modified  equinoctial  element  g  .  On  the  other  hand,  Figure  13  shows  the  evolution  over  one  orbit  of  the  out-of-plane 

modified  equinoctial  element  h .  The  nominal  value  and  the  values  due  to  errors  in  the  NKO  elements  are  shown. 
Figure  13  shows  that  the  error  on  the  out-of-plane  element  h  due  to  uncertainties  in  the  NKO  elements  is  orders  of 
magnitude  smaller  than  the  value  of  the  uncertainties.  The  same  behavior  is  noted  in  the  other  out-of-plane  modified 
equinoctial  element  k  .  Here,  the  semi-latus  rectum  shows  a  constant  error  due  to  uncertainties  in  the  NKO  elements 
of  the  same  order  of  magnitude  as  the  uncertainties  themselves.  Lastly,  the  true  longitude  exhibits  a  linearly  increasing 
error  due  to  a  0.1%  uncertainty  on  the  orbit  angular  velocity,  as  shown  in  Figure  11  for  the  RAAN  case.  After  one 
orbit,  the  absolute  error  on  the  true  longitude  is  SL  -  0.36  deg  . 
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xIO'3  xIO'3 


Figure  12.  Evolution  of  the  osculating  in-plane  modi¬ 
fied  equinoctial  element /  over  one  orbit.  Out-of- plane 
displaced  GEO  (test  case  1). 


Figure  13.  Evolution  of  the  osculating  out-of-plane 
modified  equinoctial  element  h  over  one  orbit.  Out-of- 
plane  displaced  GEO  (test  case  1). 


3.  Augmented  integrals  of  motion 

For  the  augmented  integrals  of  motion,  the  sensitivity  analysis  study  demonstrates  that  the  mapping  that  uses  these 
elements  is  robust  to  uncertainties  in  the  NKO  elements. 

For  the  sake  of  conciseness,  and  because  all  the  errors  are  constant,  the  plots  of  the  sensitivity  analysis  study  are 
not  shown  here.  Table  5  shows  a  summary  of  the  results  of  the  sensitivity  analysis  study  in  the  case  of  augmented 
integrals  of  motion. 


Table  5.  Summary  of  the  results  of  the  sensitivity  analysis  study  in  case  of  augmented  integrals  of  motion.* 


Errors  due  to 

uncertainties 

Test 

case 

Sz 

dp 

8tu 

£k 

all 

3x10 

0.2% 

0.1% 

£e 

all 

~0 

3xlO“3 

2xl0“3 

£l 

all 

0 

0 

linearly  increasing  up  to  0.2% 

Note  that  sometimes  ~0  is  used  against  0.  This  is  done  to  stress  the  difference  between  a  null  value  and  a  very  small  (numerical)  value. 
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2.1.4.  Non-Keplerian  Orbit  to  Osculating  Orbital  Elements  Map:  Summary 

A  summary  of  advantages  and  drawbacks  of  the  three  forward  mappings  investigated  in  this  study  is  shown  in 
Table  6. 

It  has  been  noted  that,  regardless  of  the  mapping  and  the  test  case  chosen,  a  spacecraft  on  a  non-Keplerian  orbit  is 
always  at  either  the  apogee  or  perigee  of  the  osculating  Keplerian  orbit.  This  can  cause  sensitivity  issues  in  case  of 
small  eccentricity.  If  the  eccentricity  of  an  orbit  if  small,  as  in  the  test  cases  presented  here,  even  small  variations  of 
the  NKO  elements  can  cause  the  eccentricity  vector  to  change  sign.  As  shown  in  Eq.  (2.23),  if  the  term 

/j  -  m 2  p2  ^Jp2  +z2  changes  sign  due  to  a  variation  in  one  of  the  NKO  elements,  the  eccentricity  vector  changes  sign. 

Similarly,  both  argument  of  perigee  and  true  anomaly  have  a  180-deg  difference  with  respect  to  their  previous  values. 
Therefore,  the  spacecraft  appears  to  instantaneously  change  its  position  in  the  osculating  Keplerian  orbit  from  apogee 
to  perigee  or  vice  versa,  without  changing  its  physical  position  in  the  Cartesian  reference  frame.  The  modified  equi¬ 
noctial  elements  are  not  affected  by  this  problem.  This  characteristic,  together  with  their  easy  formulation  and  their 
robustness  to  uncertainties  in  the  NKO  elements,  makes  the  use  of  modified  equinoctial  elements  the  best  choice  for 
the  non-Keplerian  orbit  to  osculating  orbital  elements  map. 


Table  6.  NKO  to  osculating  orbital  elements  map.  Summary  of  advantages  and  drawbacks  of  the  chosen  mappings. 


Classical  Keplerian  elements 

Modified  equinoctial 

elements 

Augmented  integrals  of 

motion 

Advantages 

•  Easy  to  understand 

•  Only  one  element  is  time  de¬ 
pendent 

•  Easy  formulation  (no  piece- 
wise  elements) 

•  Robust  to  uncertainties  in 

the  NKO  elements 

•  Robust  to  the  switch  apo¬ 
gee/perigee 

•  No  singularities 

•  Only  one  element  ( L )  is 
linearly  increasing  with 
time  if  there  are  uncertain¬ 
ties  in  m 

•  Easy  formulation 

•  No  singularities 

Drawbacks 

•  Piecewise  formulation 

•  Very  sensitive  to  uncertain¬ 
ties 

•  Issues  around  z  =  0 

•  Singularity:  Q  =  0  if  z  =  0 

•  Very  sensitive  to  the  switch 
apogee/perigee 

•  Analytical  Jacobian  unde¬ 
fined  for  3  out  of  6  elements 
because  of  the  piecewise  for¬ 
mulation 

•  Two  elements  ( Q,  co )  are  lin¬ 
early  increasing  with  time  if 
there  are  uncertainties  in  in 

•  5  elements  are  time  depend¬ 
ent 

•  Less  intuitive 

•  7  elements  instead  of  6 

•  Very  sensitive  to  the  switch 
apogee/perigee  (eccen¬ 
tricity  vector  changes  sign) 
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2.2.  Inverse  Map:  Osculating  Orbital  Elements  to  Non-Keplerian  Orbit 

After  the  forward  mappings  from  non-Keplerian  orbits  families  to  osculating  orbital  elements  described  in  Section 
2.1,  three  different  inverse  mappings  from  orbital  elements  to  non-Keplerian  orbit  geometry  have  been  generated. 
These  mappings  have  been  derived  in  closed,  analytical  form.  A  sensitivity  analysis  has  been  performed  to  understand 
the  impact  of  uncertainties  in  the  classical  orbital  elements  on  the  NKO  properties. 

The  analytical  formulation  of  the  three  inverse  mappings  is  shown  in  Section  2.2.1,  whereas  the  sensitivity  analysis 
performed  on  the  maps  is  described  in  Section  2.2.2.  Section  2.2.3  shows  the  numerical  test  cases  taken  into  account. 
Lastly,  the  three  mappings  are  compared  through  their  advantages  and  drawbacks  in  Section  2.2.4. 


2.2.1.  Mappings 

In  the  following  sub-sections,  the  inverse  mappings  from  osculating  orbital  elements  to  non-Keplerian  orbits  are 
summarized  for  the  three  sets  of  elements  considered. 


2.2.1. 1.  Classical  Keplerian  Elements 

The  map  from  osculating  Keplerian  elements  xKEP  to  Cartesian  position  and  velocity  is  here  derived.  Note  that 

the  following  map  does  not  work  in  case  the  osculating  Keplerian  orbit  is  a  parabola.  However,  no  NKOs  in  the  two- 
body  problem  are  described  by  osculating  parabolas. 

The  semi-latus  rectum  is  given  by 


p  =  a(l-e2^j 

The  radius  of  the  orbit  is 


r  =  ■ 


l  +  ecos5 


In  the  orbital  reference  frame  ,  position  and  velocity  vectors  are  given  by 


rcos$ 

r  cos  3  -  r  sin  33 

rorb  = 

rsin*9 

’  V  orb 

rsint9  +  rcos  33 

0 

0 

in  which 


r,  h  Jmp_ 

^2  2 

r  r 


dt\  l  +  e  cos  3 


ep  sin  33 
(l  +  ecos*9)2 


=  \f-esin3 


(2.106) 


(2.107) 


(2.108) 


(2.109) 
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The  rotation  matrix  R  =  R3  (—co)  — »  ri(_/)^r3(_q) 

from  orbital  reference  frame  to  Earth-Centered  Inertial 

(ECI)  is 


R  = 


cos  Q  cos  co  -  cos  i  sin  Q  sin  co 
sin  Q  cos  co  +  cos  i  cos  Q  sin  co 
sin  i  sin  co 


-  cos  Q  sin  co  -  cos  i  sin  Q  cos 

-  sin  Q  sin  ^  +  cos  i  cos  Q  cos 

sin  i  cos  co 


sin  i  sin  Q 
-  sin  i  cos  Q 
cos  i 


(2.110) 


From  Eqs.  (2.108)  and  (2.110),  the  Cartesian  position  and  velocity  vectors  in  ECI  reference  frame  are 


r  =  *r„ 

v  =  Rv„ 


(2.111) 


From  the  osculating  Keplerian  elements  to  Cartesian  position  vector  map  (Eq.  (2.1 1 1))  and  the  definition  of  the 
Cartesian  position  vector  given  in  Eq.  (2.5),  the  NKO  elements  {z,p\  are  derived  as  follows. 


Z  =  K 


i  +  e  cos  1 9 


(sin  i  sin  co  cos  3  +  sin  i  cos  &>sin  3)  =  - 


+  e  cos  3 


sin  i  sin 


(co  +  3) 


(2.112) 


a(l-e2^j 
1  +  ecos  3 


yj(Rn  cos  3  +  Rn  sint9)2  +(R21  cost9  +  /?22  sin$)2 


(2.113) 


After  some  algebraic  manipulations  of  Eq.  (2.1 13),  the  formulation  for  the  orbital  radius  p  is 

a(  l-e2)  i - - - - — 

/7=2(l  +  ecos^)^3  +  COS^  +  2C°S^2^  +  l9^Si 


sin2  i 


(2.114) 


For  what  concerns  the  orbital  angular  velocity  m  ,  let  us  compute  first  the  magnitude  of  the  velocity  vector.  Be¬ 
cause  the  rotation  matrix  R  shown  in  Eq.  (2.1 10)  does  not  change  the  magnitude  of  the  vector,  it  is  easier  to  compute 
the  magnitude  of  the  velocity  vector  in  the  orbital  reference  frame  (Eq.  (2.108)). 


v  =  ^rcosi9-rsini9i9^  +(rsin$  +  r  cos  33^  =  ...  =  4 f 2  +r232 
Using  Eq.  (2.109)  and  after  some  algebraic  manipulation,  Eq.  (2.115)  becomes 


v  =  J~(  1  +  +  2e  cos  3) 


(2.115) 


(2.116) 


Therefore,  the  formulation  of  the  orbital  angular  velocity  zu  is  given  by 

^  -2(l  +  £cos$) 


P  ]la3(l-e2f 


l  +  e  +2^  cos  3 


1 3  +  cos(2/)  +  2cos(2(&>  +  -9)) 


sin2  i 


(2.117) 
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Summary 

The  inverse  map  from  osculating  classical  Keplerian  elements  to  NKOs  elements  described  in  Section  2.2. 1.1  is 
summarized  as 

a(l-e2) 

z  =  — - —  sin  i  sin  (co  + 

l  +  ecosi9 

i - 

2(l  +  ecosi9)  * 

I  Li  I  1  +  e2  +  2ecosi9 

- — - v2(l  +  ecos*9)  - — — - r- - tt — — 

v  a3  (l-^2)  3  + cos  (2/) +  2  cos  (2  (&>  +  $))  sin2  i 


$  +  cos 


(2/)  +  2cos(2(&>  +  si 


sin2  i 


(2.118) 


2.2. 1.2.  Modified  Equinoctial  Elements 

The  map  from  osculating  modified  equinoctial  elements  xMEE  to  Cartesian  position  and  velocity  is  now  derived. 
First,  let  us  define  the  following  auxiliary  variables,  which  will  be  useful  in  the  mathematical  derivation  to  shorten 
the  equations. 


q  =  1  +  /  cos  L  +  g  sin  L 
r  =  p/q 
s  —  1  +  h  +  k 
J3  =  h2-k2 


(2.119) 


From  Broucke  and  Cefola,23  and  after  some  algebraic  manipulation  due  to  the  difference  in  the  elements  consid¬ 
ered,  the  formulation  of  the  Cartesian  position  vector  is  given  as  from  Betts  (Eq.  6.42,  page  266). 24 


r 

r  =  - 


cos  L  +  p  cos  L  +  2  hk  sin  L 
sin  L- P sin  L  +  2 hk  cos  L 
2(/zsinL-kcosL) 


(2.120) 


From  the  formulation  of  the  position  vector  in  Eq.  (2.120),  the  Cartesian  velocity  vector  is  given  by  Eq.  (2.121), 
in  which  dL/dt  -  y[jup(q/ p )2 ,  as  from  De  Pascale  and  Vasile.25 


v  = 


dr  dL 
8L  dt 


(2.121) 


Therefore,  the  formulation  of  the  Cartesian  velocity  vector  given  by  Eq.  (2.121)  is  the  same  as  from  Betts  (Eq. 
6.43,  page  266). 24 


v  =  i> 


SMp 


-  sin  L  -  P  sin  L  +  2 hk  cos  L  -  g  +  2  fhk  -J3g 
cos  L  -  p  cos  L  -  2 hk  sin  L  +  /  -  2 ghk  -J3f 
2[h  cos  L  +  k  sin  L  +  fh  +  gk) 


(2.122) 
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From  the  osculating  modified  equinoctial  elements  to  Cartesian  state  map  (Eqs.  (2.120)  and  (2.122))  and  the  defi¬ 
nition  of  the  Cartesian  position  vector  given  in  Eq.  (2.5),  the  NKO  elements  [z,p]  are  derived  as  follows. 


2r 


Z  =  rz  =  — (/zsinL-&cosL) 

P  =  frx  +ry  =  -^|(cosL  +  /?cosL  +  2MsinL)2  +(sinL-/?sinL  +  2McosL)2 


(2.123) 

(2.124) 


After  some  algebraic  manipulations  of  Eq.  (2.124),  the  formulation  for  the  orbital  radius  p  is 


P  = 


-yjl  +  (h2  +fc2)2  +2y9cos(2L)  +  4Msin(2L) 


(2.125) 


For  the  orbital  angular  velocity  m  ,  let  us  compute  first  the  magnitude  of  the  velocity  vector.  From  Eq.  (2.122), 
the  magnitude  of  the  velocity  vector  is 


v  —  ,  j- —  +  f  +  g  +2  f  cos  L  +  2g  sin  L 


(2.126) 


Therefore,  the  formulation  of  the  orbital  angular  velocity  m  is 


vj  -  - 


s  I p  I  1  +  /  +2y cos L  +  2g sin L 

P  V 1  +  ( h2  +  k2  )2  +  2 J3 cos (2L)  +  4 hk  sin  (2 L) 


(2.127) 


Summary 

The  inverse  map  from  osculating  modified  equinoctial  elements  to  NKOs  elements  described  in  Section  2.2. 1.2  is 
summarized  as 


in  which 


2r 


z  =  —  (hsinL  —  k  cos  L ) 


» =  -  ^1  +  (V  +  T2  f  +  2/?  cos  (2L)  +  4hk  sin  (2L) 


s  I p  I  1  +  f  +  g  +2 f  cos  L  +  2g  sin 
r  Ip  ]jl  +  (h2  +  k2 )2  +  2p cos (2L)  +  4hk sin (2 L) 


1  +  /  cos  L  +  g  sin  L 
<  s  =  l  +  h2  +k2 


(2.128) 


(2.129) 


J3  =  h2  -k2 
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2.2.I.3.  Augmented  Integrals  of  Motion 


The  map  from  osculating  augmented  integrals  of  motion  xAIoM  to  Cartesian  position  and  velocity  is  now  derived. 

Let  us  recall  here  the  following  expressions  that  link  the  integrals  of  motion  with  the  Cartesian  position  and  velocity, 
which  will  be  useful  for  the  mathematical  derivation  of  the  inverse  mapping. 


h  =r  xv 


(2.130) 


rJ 


(2.131) 


The  magnitude  of  the  orbital  angular  momentum  is 

1*1  =  \[pM  = 

Therefore,  the  semi -major  axis  can  be  retrieved  by  Eq.  (2.132). 


h 


a  =  - 


M) 


Reminding  us  that  the  instantaneous  position  vector  along  the  NKO  corresponds  always 
perigee  of  the  osculating  Keplerian  orbit,  the  magnitude  of  the  position  vector  is  given 
plus/minus  sign  is  referred  to  the  case  of  apogee/perigee,  respectively. 


*(!  +  ,)  = 


(l±e) 


(2.132) 


(2.133) 

with  either  the  apogee  or  the 
by  Eq.  (2.134),  in  which  the 


(2.134) 


The  direction  of  the  position  vector  is  the  same  as  that  of  the  eccentricity,  as  shown  in  Eq.  (2.131).  However,  there 

/ 

is  an  uncertainty  on  the  sign  of  the  position  unit  vector  which  is  related  to  sgn 


v  o 


Kt*  r) 


r  =  < 


—e  if 

e  if 


v_r 

v/'  n 
V  O 


<  0  (apogee) 
>  0  (perigee) 


r) 

From  Eqs.  (2.134)  -  (2.135),  the  formulation  of  the  Cartesian  position  vector  is  given  by 


r  =  < 


/*(!-«) 


l+e) 


e  (apogee) 


e  (perigee) 


(2.135) 


(2.136) 


The  ambiguity  on  the  sign  in  Eq.  (2.136)  due  to  the  uncertainty  on  the  apogee  or  perigee  position  can  be  solved  by 
using  the  information  given  by  the  true  longitude.  That  is,  the  difference  between  the  true  longitude  and  the  longitude 
of  perigee  ( co  +  Q )  can  be  either  null  (if  the  spacecraft  is  at  the  perigee  of  the  osculating  orbit)  or  1 80  deg  (if  the 
spacecraft  if  at  the  apogee  of  the  osculating  orbit).  The  longitude  of  perigee  is  given  by  the  eccentricity  unit  vector  as 
follows. 
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co  +  Q  =  arctan 


'e  e' 


=  arctan  (ex,ey) 


(2.137) 


Therefore,  apogee  and  perigee  are  defined  as  follows. 


J apogee  if  mod  (a>  +  Q  +  n,  2 71)  -  L 
[perigee  if  co  +  Q  =  L 


(2.138) 


in  which  all  the  angles  are  considered  to  lie  in  the  open  interval  [0,  2tt)  . 

From  the  definition  in  Eq.  (2.138),  the  formulation  of  the  Cartesian  position  vector  in  Eq.  (2.136)  can  be  rewritten  as 


r  = 


\m  „  /  t  \  \ 

— — -  e  if  mod  arctan  ( e  ,  e  )  +  n,  2n  =  L 

ju(l-e )  v  v  7  ’ 

- e  if  arctan  (ex ,  ey  )  =  L 


(2.139) 


From  Eq.  (2.130),  because  the  radial  component  of  the  velocity  is  always  zero  for  the  NKOs  considered,  the  mag¬ 
nitude  of  the  velocity  vector  is  simply  given  by 


r 


P^j(l-e)  if  mod^arctan[^,^)  +  ^-,2^-j  =  L 
TjTjr(l  +  e)  if  arctan (ex,ey)  =  L 


(2.140) 


The  direction  of  the  velocity  is 


v =  hxf  =  < 


-  [h  x  e  j  if  mod  ^arctan  [ex  ,ey^  +  7r,  2n ^  =  L 
[h  x  e  j  if  arctan  (ex ,  ey )  =  L 

Therefore,  the  expression  of  the  velocity  vector  is  given  by 

-|Mi(l-^)(^X^)  ^  mod^arctan^,^  )  +  ;r,2;r)  =  L 
Irll 

jj^jj(l  +  e){hxe^  if  arctan [ ex, ey  )  =  L 


(2.141) 


v  =  < 


(2.142) 


From  the  osculating  augmented  integrals  of  motion  to  Cartesian  position  vector  map  (Eq.  (2.139)),  and  the  defini¬ 
tion  of  the  Cartesian  position  vector  given  in  Eq.  (2.5),  the  parameters  z  and  p  are  derived  as  follows. 


z  =  r  = 


, 

- — -ez  if  mod^arctan[^x,^)  +  ^-,2^-^  =  L 


(2.143) 


p(\  +  e) 


ez  if  arctan  (ex ,  ey )  =  L 
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=i 


/f(i_g)^gj/gl2+(gy/e)2  ifmod(arctan(^’e,)+;r’2;r)=i 

j^ii2  I - - - — 

77~^yl(ex/e)2  +{ey/e)  if  arctan (ex,ey)  =  L 


(2.144) 


//(l  +  e) 


From  the  osculating  augmented  integrals  of  motion  to  Cartesian  velocity  vector  magnitude  map  (Eq.  (2.140)),  and 
the  definition  of  the  Cartesian  velocity  vector  given  in  Eq.  (2.6),  the  parameter  ur  is  derived  as  follows. 


vj  -  —  -  — -  -  \ 
P  rp 


i 


ll*H 

A 2  ( \+ef 


^(ex/e)2+(ey/ef 


if  mod  (arctan  (ex ,  ey )  +  n,  2n  )  =  L 


if  arctan  =  L 


(2.145) 


Summary 

The  inverse  map  from  osculating  augmented  integrals  of  motion  to  NKOs  elements  described  in  Section  2. 2. 1.3  is 
summarized  as 


z  =  < 


p\-e) 


if  mod  (arc tan  (ex ,  ey )  +  n,  2jt^  —  L 
if  arctan  (er,cy)  =  L 


I  =  L 


in  =  < 


p\+e) 

■■  ii2  | - — 

- — -^|(ex/e)2  +(ey/e)  ^  mod(arctan(ex,ey)  +  ;r,2;r)  = 

^7jef7{7Jef  if  arctan  (ex,ey)  =  L 

if  mod  (arctan  (ex ,  ey )  +  n,  2n  j  =  L 


(2.146) 


/j(l  +  e) 
^(1-c) 


^2(l  +  e)2 


^j4+[ejef 

_ 1 _ 

^Jef+{ejef 


if  arctan  [ex ,  ey )  =  L 


A.  Pelonit5|^femp|(^Ki^proved  for  pub|jc  re|ease:  distribution  unlimited.  Pa§e  39 


Osculating  Keplerian  Elements  for  Highly  Non-Keplerian  Orbits  -  Final  technical  report 


2.2.2.  Sensitivity  Analysis 

A  sensitivity  analysis  has  been  performed  in  order  to  study  the  behavior  of  the  mappings  in  the  presence  of  errors 
in  the  osculating  orbital  elements.  Only  a  numerical  sensitivity  analysis  has  been  carried  out  by  computing  the  values 
of  the  osculating  elements  from  the  perturbed  input  elements  and  then  computing  the  errors  with  respect  to  the  nominal 
values. 

Denoting  with  y  the  set  of  elements  chosen  to  describe  the  osculating  Keplerian  orbit,  the  parameters  used  to  test 
the  sensitivity  of  the  mapping  to  uncertainties  in  the  osculating  orbital  element  y.  are  as  follows. 

•  Absolute  error  of  the  vertical  displacement  due  to  uncertainties  in  the  input  elements.  The  absolute  error 
has  been  chosen  against  the  relative  error  because  all  the  test  cases  under  study  are  characterized  by  small 
vertical  displacements. 


sz^  :=\z-z\  =  \z-(z  +  Sz)\  =  \Sz\  (2.147) 

•  Relative  error  of  the  orbit  radius  due  to  uncertainties  in  the  input  elements.  The  relative  error  has  been 
chosen  against  the  absolute  error  because  all  the  test  cases  under  study  are  characterized  by  large  orbit 
radii. 


\p-p\_  \p-(p  +  sp)\_  \Sp\ 

£p,yt  ■_  ~  ~ 

P  P  P 

Relative  error  of  the  orbital  angular  velocity  due  to  uncertainties  in  the  input  elements. 


(2.148) 


\m-m\  m 


£m,yi  ■ 


-(zn  +  Szn)\  |  St 


m\ 

m 


(2.149) 


2.2.3.  Numerical  Test  Cases 

The  same  two  sets  of  test  cases  described  in  Section  2.1.3  have  been  considered.  Therefore,  the  description  of  such 
test  cases  is  not  repeated.  Section  2.2.3. 1  shows  a  validation  of  the  inverse  maps,  whereas  is  Section  2. 2. 3. 2  the  sen¬ 
sitivity  analysis  on  the  inverse  maps  is  discussed. 


2.2.3. 1.  Forward  and  Inverse  Mappings 

All  the  results  of  the  test  cases  described  in  Section  2. 1.3.1  have  been  numerically  verified  against  the  nominal 
values.  That  is,  the  osculating  orbital  elements  have  been  obtained  through  the  forward  mappings  starting  from  the 

nominal  NKO  elements  [z,p,mf  .  Then,  the  inverse  mappings  are  used  to  obtain  the  original  NKO  elements.  The 

errors  between  the  nominal  values  of  the  NKO  elements  and  those  obtained  after  the  conversion  give  an  estimate  of 
the  accuracy  of  both  forward  and  inverse  mappings. 

Figures  14,  15,  and  16  show  the  absolute  errors  after  forward  and  inverse  mappings  in  out-of-plane  displacement, 
orbit  radius,  and  orbital  angular  velocity,  respectively.  It  is  clear  that  the  errors  after  the  forward  and  inverse  mappings 
are  negligible.  This  validates  both  the  forward  and  inverse  mappings. 
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Figure  14.  Absolute  error  in  the  out-of-plane  dis¬ 
placement  after  forward  and  inverse  mappings. 


Figure  15.  Absolute  error  in  the  orbit  radius  after 
forward  and  inverse  mappings. 
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Figure  16.  Absolute  error  in  the  orbital  angular  velocity  after  forward  and  inverse  mappings. 
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2.2.3.2.  Sensitivity  Analysis 

The  same  test  cases  taken  into  account  for  the  sensitivity  analysis  of  the  forward  mappings  as  described  in  Section 
2. 1.3.2  have  been  considered  here. 

The  numerical  sensitivity  analysis  has  been  performed  as  described  in  Section  2.2.2. 

The  errors  of  the  osculating  classical  Keplerian  elements  are  considered  fixed  and  constant,  as  follows. 

^=10-3 
st  =  10-3 

£•  =  1  tr3  ■  7T 

£n  =10”3  • 7Z 
£a  =10^-71 
£g  =  10“3  -7Z 

The  errors  of  the  osculating  modified  equinoctial  elements  are  considered  fixed  and  constant,  as  follows. 


(0.1%) 


(0.1%  of  the  maximum  possible  error  -  0.1 8 deg) 
(0.1%  of  the  maximum  possible  error  -  0.1 8 deg) 
(0.1%  of  the  maximum  possible  error  -  0.1 8 deg) 
(0.1%  of  the  maximum  possible  error  -  0.1 8  deg) 


(2.150) 
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(0.1%) 


(0.1%  of  the  maximum  possible  error  -  0.1 8  deg) 


(2.151) 


No  sensitivity  analysis  has  been  performed  in  the  case  of  augmented  integrals  of  motion.  This  is  because,  as  shown 
in  Eq.  (2.146),  the  augmented  integrals  of  motion  to  NKO  elements  map  is  extremely  sensitive  to  errors  in  both  ec¬ 
centricity  vector  and  true  longitude.  The  issue  is  not  only  the  piecewise  formulation  per  se  but  the  fact  that  the  formu¬ 
lation  takes  into  account  only  two  points  in  the  orbit,  which  are  the  apogee  and  perigee.  If  an  error  in  either  eccentricity 
vector  or  true  longitude  moves  the  point  from  the  apogee/perigee,  the  inverse  map  fails. 


Results  and  discussion 

The  sensitivity  analysis  study  demonstrates  that  the  planar  elements  (i.e.  )  of  both  mappings  are  robust  to 

uncertainties  in  the  osculating  orbital  elements.  Both  mappings  of  the  out-of-plane  displacement  z  are  robust  to  un¬ 
certainties  in  the  in-plane  osculating  elements.  However,  the  out-of-plane  displacement  is  very  sensitive  to  errors  in 
the  out-of-plane  osculating  elements  (i.e.  inclination,  if  KEP  are  used,  or  out-of-plane  MEE  h,k  ).  This  is  understand¬ 
able  if  the  Keplerian  orbit  has  a  small  inclination  and  a  large  semi -major  axis  ( a  ~  rGE0  in  this  case). 

For  the  sake  of  conciseness,  only  the  most  interesting  plots  of  the  sensitivity  analysis  study  are  shown  here.  Table 
7  and  Table  8  show  a  summary  of  the  numerical  results  of  the  sensitivity  analysis  study  in  the  case  of  KEP  and  MEE, 
respectively.  The  three  test  cases  taken  into  account  are  those  described  in  Eqs.  (2.99)  -  (2.101). 
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Table  7.  Summary  of  the  results  of  the  sensitivity  analysis  study  in  the  case  of  classical  Keplerian  elements.* 


Uncertainty 

Test  case 

Sz 

£P 

1 

~  150  m  (0.1%) 

5a 

2 

0 

0.1% 

0.15% 

3 

<0.1  m  (<0.1%) 

1 

~  150  m  (0.1%) 

0.1% 

~  0.2% 

5e 

2 

0 

0.101% 

The  three  test  cases  are 

3 

<0.1  m  (<0.1%) 

0.1% 

numerically  different. 

1 

~  130  km 

1.6x10“3% 

1.6xl0“3% 

Si 

2 

<130  km 

<5x10^% 

<5x10^% 

3 

~  130  km 

5x10^% 

5x10^% 

1 

~0 

0 

0 

sn 

2 

0 

~0 

~0 

3 

~0  (0.001%) 

0 

0 

1 

~  0.7  m  (5x10^%) 

~0 

~0 

5(0 

2 

0 

The  three  test  cases  are 

The  three  test  cases  are 

3 

~  10‘3  m  (0.002%) 

numerically  different. 

numerically  different. 

1 

~  0.7  m  (5x10”4%) 

0 

~0 

59 

2 

0 

5xl0“6% 

10“5% 

3 

~  10‘3  m  (0.002%) 

0 

0 

Note  that  sometimes  ~0  is  used  against  0.  This  is  done  to  stress  the  difference  between  a  null  value  and  a  very  small  (numerical)  value. 


Figure  17  shows  the  evolution  of  the  out-of-plane  displacement  z  over  one  orbit  obtained  from  the  classical  Kep¬ 
lerian  elements.  The  three  test  cases  are  shown.  Both  the  nominal  value  and  a  0.1%  error  in  the  inclination  are  consid¬ 
ered.  It  is  shown  how  a  small  error  in  inclination  can  cause  a  large  error  in  the  out-of-plane  displacement.  Figure  18 
shows  the  relative  error  evolution  of  the  orbit  radius  p  over  one  orbit  obtained  from  the  classical  Keplerian  elements. 

The  three  test  cases  are  shown.  A  0.1%  error  in  the  inclination  is  considered.  It  is  shown  how  an  error  in  inclination 
does  not  affect  the  orbit  radius. 


A.  Pelonit5|^femp|(^Ki^proved  for  pub|jc  re|ease:  distribution  unlimited.  Pa§e  43 


Osculating  Keplerian  Elements  for  Highly  Non-Keplerian  Orbits  -  Final  technical  report 


x1(T 


16uGOO®GOG0GGG0GGO&GGGeGGO000000000G©000000000O0©0®' 


-©-Type  1  NKO 
Type  3  NKO 

a-  Type  1  NKO  (very-small  displacement) 


12 


OO 


t  8 
<D 


0.2 


0.4  0.6 

Orbit  fraction 


0.8 


Figure  17.  Out-of-plane  displacement.  Nominal  case 
and  0.1%  error  in  inclination. 


Figure  18.  Relative  error  of  the  orbit  radius  in  case  of 
0.1%  error  in  inclination. 


Table  8.  Summary  of  the  results  of  the  sensitivity  analysis  study  in  the  case  of  modified  equinoctial  elements.* 


Uncertainty 

Test  case 

1 

~  150  m  (0.1%) 

Sp 

2 

0 

0.1% 

0.15% 

3 

<0.1  m  (<0.1%) 

1 

~  150  m  (0.1%) 

Sf 

2 

0 

<0.1% 

<  0.2% 

3 

<0.1  m  (<0.1%) 

1 

~  150  m  (0.1%) 

<5g 

2 

0 

<0.1% 

<  0.2% 

3 

<0.1  m  (<0.1%) 

Sh 

1 

<85  km 

<  10  3% 

<  10  3% 

2-3 

<2x10^% 

<2x10^% 

8k 

1 

<85  km 

A 

© 

U) 

A 

© 

U) 

2-3 

<2x10^% 

<2x10^% 

1 

-  0.7  m  (5xK<%) 

~0 

~0 

8L 

2 

0 

5xl0“6% 

10“5% 

3 

~0  (8x10_5%) 

-0 

-0 

Note  that  sometimes  ~0  is  used  against  0.  This  is  done  to  stress  the  difference  between  a  null  value  and  a  very  small  (numerical)  value. 
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Figure  19  shows  the  relative  error  evolution  of  the  orbit  radius  p  over  one  orbit  obtained  from  the  MEE  with  a 
0.001  error  in  the  in-plane  element  /  .  The  three  test  cases  are  shown.  It  is  shown  that  the  relative  error  on  the  orbit 

radius  is  at  most  the  same  order  of  magnitude  as  the  initial  uncertainty  itself.  Figure  20  shows  the  evolution  of  the  out- 
of-plane  displacement  z  over  one  orbit  obtained  from  the  MEE.  The  three  test  cases  are  shown.  Both  the  nominal 
value  and  a  0.001  error  in  the  out-of-plane  element  h  are  considered.  It  is  shown  how  a  small  error  in  the  out-of-plane 
element  can  cause  a  large  error  in  the  out-of-plane  displacement. 


0  0.2  0.4  0.6  0.8  1 

Orbit  fraction 


0  0.2  0.4  0.6  0.8  1 
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Figure  19.  Relative  error  in  orbit  radius  in  case  of 
0.001  error  in  the  in-plane  element  /. 


Figure  20.  Out-of-plane  displacement.  Nominal  case 
and  0.001  error  in  the  out-of- plane  element  h . 


2.2.4.  Osculating  Orbital  Elements  to  Non-Keplerian  Orbit  Map:  Summary 

A  summary  of  advantages  and  drawbacks  of  the  inverse  mappings  investigated  in  this  study  is  shown  in  Table  9. 


Table  9.  Osculating  orbital  elements  to  non-Keplerian  orbit  map.  Summary  of  advantages  and  drawbacks  of  the  chosen 
mappings. 


Classical  Keplerian  elements 

Modified  equinoctial 

elements 

Augmented  integrals  of 

motion 

Advantages 

•  Easy  to  understand 

•  Easy  formulation 

•  Robust  to  uncertainties 

•  Easy  formulation 

•  Robust  to  uncertainties  on 
in-plane  elements 

Drawbacks 

•  Out-of-plane  displacement 
very  sensitive  to  uncertainties 
on  inclination 

•  Out-of-plane  displacement 
very  sensitive  to  uncertain¬ 
ties  on  out-of-plane  ele¬ 
ments 

•  Piecewise  formulation 

•  Extremely  sensitive  to  er¬ 
rors  in  both  eccentricity  vec¬ 
tor  and  true  longitude  to  the 
point  that  the  map  fails  if 
there  is  an  error  in  either  ec¬ 
centricity  or  true  longitude 
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2.3.  Determination  of  Spacecraft  Thrust-Induced  Acceleration 

The  magnitude  and  direction  of  the  spacecraft  thrust-induced  acceleration  used  to  generate  families  of  non -Kep¬ 
lerian  orbits  have  been  extracted  from  the  inverse  mapping  from  the  orbital  elements  discussed  in  Section  2.2.  This 
has  been  achieved  by  using  the  NKO  geometric  properties  derived  from  the  inverse  mapping  as  inputs  to  existing 
analysis  of  the  dynamics  of  non-Keplerian  orbits.2  This  analysis  links  the  thrust-induced  acceleration  to  the  non-Kep- 
lerian  orbit  geometry. 

The  following  sub-sections  show  the  analytical  expressions  obtained  for  the  thrust-induced  acceleration  if  either 
the  classical  Keplerian  elements  or  the  modified  equinoctial  elements  are  used  to  describe  the  osculating  Keplerian 
orbit.  Because  of  the  very  poor  robustness  of  the  augmented  integrals  of  motion  shown  in  Section  2.2. 1.3  for  the 
inverse  mapping,  these  elements  have  not  been  considered  in  the  study  of  the  thrust-induced  acceleration. 

2.3.1.  Classical  Keplerian  elements 

From  Eq.  (2.1)  and  using  the  mappings  shown  in  Eq.  (2.1 18),  the  thrust-induced  acceleration  and  the  thrust  pitch 
angle  are  obtained  as 


\a 


(2.152) 


(2.153) 


in  which 


(2.154) 


2.3.2.  Modified  equinoctial  elements 


From  Eq.  (2.1)  and  using  the  map  shown  in  Eq.  (2.128),  the  thrust-induced  acceleration  and  the  thrust  pitch  angle 


are  obtained  as 


(2.156) 


in  which  q  ,  s  ,  and  /?  are  defined  in  Eq.  (2.119),  whereas  <r2  is 


cr2  =  \  +  {h2  +k2)  +2/?cos(2L)  +  4/zksin(2L) 


(2.157) 
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2.3.3.  Numerical  Test  Cases 


The  same  two  sets  of  test  cases  described  in  Section  2.1.3  have  been  considered  here.  Therefore,  the  description 
of  such  test  cases  is  not  repeated. 

As  noted  earlier,  the  NKOs  generated  by  the  thrust-induced  acceleration  shown  in  Eq.  (2.1)  correspond  to  circular 
orbits  displaced  above  the  central  body,  viewed  from  an  inertial  frame  of  reference.  Therefore,  a  displacement  below 
the  central  body  is  not  explicitly  taken  into  account.  This  is  because  the  tangent  of  the  thrust  pitch  angle,  rather  than 
the  pitch  angle  itself,  is  given.  Therefore,  there  is  an  ambiguity  in  the  thrust  pitch  angle.  This  can  be  solved  algorith¬ 
mically  by  checking  the  sign  of  the  out-of-plane  displacement,  so  that 


ftan  1  (tan  or)  ifz>0 
[tan-1  (tan  a)  +  7r  if  z  <  0 


(2.158) 


The  acceleration,  in  terms  of  magnitude  and  pitch  angle,  resulting  from  both  Eqs.  (2. 152)  -  (2. 153)  and  Eqs.  (2. 155) 
-  (2.156),  has  been  validated  using  two  separate  checks.  In  the  first  check,  the  values  of  acceleration  magnitude  and 
pitch  angle  are  compared  with  those  given  directly  by  using  Eq.  (2.1).  The  values  of  both  acceleration  magnitude  and 
the  direction  of  the  acceleration  have  been  verified  for  all  test  cases  considered  with  both  sets  of  elements.  The  errors 
in  the  magnitude  and  direction  of  the  acceleration  are  ||<?a||  <  2xl0-9  mm/s2  and  \Sa\  <  10"10  rad,  respectively.  The 

second  check  has  been  carried  out  by  propagating  the  equations  of  motion  with  the  propulsive  acceleration  given  by 
Eqs.  (2.152)  -  (2.153)  and  Eqs.  (2.155)  -  (2.156),  respectively.  The  relative  errors  in  position  vector  and  velocity 
vector  between  the  initial  state  and  the  state  after  one  orbital  period  are  used  to  validate  the  acceleration.  A  variable 
order  Adam-Bashford-Moulton  PECE  solver  (as  implemented  in  MATLAB  ode  113)  with  relative  and  absolute  toler¬ 
ances  equal  to  10-10  and  10-12 ,  respectively,  has  been  used  for  the  propagation  of  the  equations  of  motion.  The  errors 
in  the  final  position  and  velocity  are  ||<5r||/||r||  <  7xl0-8  and  ||£v||/||v||  <  10-7 ,  respectively. 


Figure  21  shows  the  Type  1  NKOs  described  in  Eq.  (2.93)  together  with  the  acceleration  vector  required.  The 
acceleration  is  clearly  shown  to  be  always  directed  along  the  vertical  axis.  Figure  22  shows  the  in-plane  Type  3  NKOs 
described  in  Eq.  (2.97)  together  with  the  acceleration  vector  required.  In  this  case,  it  is  clear  how  the  propulsive 
acceleration  is  used  to  balance  the  gravitational  attraction  of  the  primary  body.  In  the  case  of  p<  rGE0  ,  for  instance, 

the  NKO  is  slower  than  a  Keplerian  orbit  with  the  same  radius.  The  spacecraft,  in  this  case,  must  thrust  in  the  radial 
direction  to  balance  the  gravitational  force  and  be  able  to  move  at  the  desired  angular  velocity. 
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x  [km]  xio4 


Figure  21.  Type  1  NKOs  and  acceleration  (not  to  Figure  22.  In-plane  displaced  Type  3  NKOs  and  ac- 

scale)  needed  with  GEO  shown  (dashed  cyan  line).  celeration  (not  to  scale)  needed  with  GEO  shown  (dashed 

cyan  line). 

As  noted  in  Section  2.1.3. 1  when  the  test  cases  used  to  verify  the  mappings  have  been  introduced,  these  cases  have 
been  chosen  for  ease  of  illustration  with  the  only  purpose  to  verify  the  validity  of  the  mappings  and  study  the  behavior 
of  the  osculating  elements.  That  is,  such  test  cases  were  not  meant  to  be  examples  of  feasible  missions.  This  is  evident 
if  the  magnitude  of  acceleration  required  is  considered.  The  accelerations  required  for  the  test  cases  shown  in  Section 
2. 1.3.1  have  magnitudes  ||a||  g  [60,200]  mm/s2  ,  which  clearly  are  too  large  for  a  conventional  low-thrust  spacecraft. 

However,  the  test  cases  described  in  Section  2. 1.3.2  for  the  study  of  the  sensitivity  analysis  have  been  chosen  as 
feasible  mission  scenarios.  The  accelerations  required  for  this  second  set  of  test  cases  have  magnitudes 

||« ||  g  [4  x  1(T4 ,  2.3]  mm/ s2  .  Considering  a  1 000-kg  spacecraft,  such  acceleration  can  be  easily  converted  into  a  max¬ 
imum  thrust  7^  g[4x10~\  2.3 J  N  that  are  feasible  accelerations  for  near-term  low-thrust  technology.26 
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3.  Generalized  Model:  No  Assumptions  on  Orbital  Plane 

The  mappings  derived  in  Section  2  are  related  to  non-Keplerian  orbits  whose  orbital  plane  is  parallel  to  the  X  —Y 
plane  in  a  Cartesian  reference  frame  (i.e.  the  equatorial  plane,  in  this  work).  A  generalized  mapping  for  the  direct 
problem  is  here  investigated,  for  example  where  the  orbit  plane  can  be  rotated  relative  to  the  Earth’s  equatorial  plane. 


Figure  23.  Sketch  of  a  displaced  orbit  with  thrust-induced  acceleration.  Generalized  model. 


Two  extra  elements  are  introduced  for  the  description  of  the  non-Keplerian  orbit  in  the  three-dimensional  space, 
as  shown  in  Figure  23.  These  are 

f  j  Inclination  of  the  non-Keplerian  orbit 

[g  Right  ascension  of  the  ascending  node  of  the  non-Keplerian  orbit 

Note  that  the  symbols  {  j,  £}  have  been  used  for  the  description  of  the  NKO  elements  inclination  and  right  ascen¬ 
sion  of  the  ascending  node  (RAAN)  because  these  elements  do  not  have  to  be  confused  with  the  osculating  inclination 
and  RAAN,  respectively.  Therefore,  the  set  of  generalized  NKO  elements  is  xNKO  =  [z,/?,ezj,  . 


3.1.  Forward  Map:  Non-Keplerian  Orbit  to  Osculating  Orbital  Elements 


The  formulation  of  the  Cartesian  position  and  velocity  vectors  shown  in  Eqs.  (2.5)  -  (2.6)  can  be  generalized, 
considering  inclination  and  RAAN,  as 


>o(cos4:cos(s7^)-cos  jsin^sin(67t))  +  zsin  jsin^ 
yo(sin^cos(^)  +  cos  j  cos  ^  sin  zsin  j  cos^ 

p  sin  jsin(s7t)  +  zcos  j 


(3.2) 


v  =  mp 


-cos^sin(sjr)-cos  j  sin  cos  (W) 
-sin4:sin(^)  +  cos  j  cos^cosfW) 
sin  jcos(ftTt) 


(3.3) 
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The  orbital  angular  momentum  h  can  be  computed  by  using  the  aforementioned  expressions  of  position  and  ve¬ 
locity  as 


h-rxv-mp 


p  sin  jsin^-zcos^cosfW/^  +  zcos  j  sin^sin^mt) 
-z  sin  cos  (W)-/?  sin  jcos^-zcos  7  cos  sin  (W) 
p  cos  j-zsin  jsin(mt) 


(3.4) 


Following  the  same  procedure  used  to  obtain  the  mappings  in  the  simplified  case,  the  generalized  mappings  from 
NKO  elements  to  osculating  orbital  elements  are  obtained.  Because  of  the  very  poor  robustness  of  the  osculating 
AIoM  for  the  inverse  mapping,  these  elements  have  not  been  considered  in  this  study. 

In  the  following  sub -sections,  the  generalized  forward  mappings  from  non -Keplerian  orbit  to  osculating  orbital 
elements  are  derived  for  the  two  sets  of  elements  considered. 


3.1.1.  Classical  Keplerian  Elements 

Starting  from  the  Cartesian  state  vector  (Eqs.  (3.2)  -  (3.3)),  the  mapping  from  generalized  NKO  elements 
\z,p,vj,j,E^  t0  osculating  classical  Keplerian  elements  is  obtained  as  described  in  Algorithm  4.2  from  Reference 
22. 

7.  Inclination  ( i  ) 


1  =  arccos 


(ft.z)= 


arccos 


p  cos  j  -  z  sin  7  sin  (art) 


t 


p2  +z2 


(3.5) 


Note  that,  if  the  NKO  orbital  plane  is  parallel  to  the  equator,  Eq.  (3.5)  becomes  the  following,  which  is  the  de¬ 
scription  already  obtained  in  the  simplified  case  (Eq.  (2.23)). 


z'h  =  arccos  {h  ( j  -  0)  •  Z  j  =  arccos 


f  \ 

P 


1 


p2  +  z2 


(3.6) 


2.  Right  Ascension  of  the  Ascending  Node  (Q  ) 

The  line  of  the  nodes  shall  be  defined  in  order  to  find  the  expression  for  the  right  ascension  of  the  ascending  node. 
The  line  of  the  nodes  is  defined  as 


N. 

Z  xh 


p  sin  y  cos<f  +  z(cos  j  cos^sin^mt)  +  sin4:cos(W)) 
p  sin  j  sin  <^  +  z  (cos  j  sin  £  sin  (s?t)  -cos^cos(G7t)) 
0 


(3.7) 


in  which 

N  =  ^z2  cos2  (^)  +  (/?sin  j  +  zcos  jsin(s7r))2  (3.8) 

It  is  important  to  note  that  the  line  of  the  nodes  is  not  defined  in  the  case  of  a  planar  orbit  with  the  orbital  plane 
parallel  to  the  equatorial  plane  (i.e.  z  =  j  =  0).  Therefore,  the  line  of  the  nodes,  in  this  case,  is  arbitrarily  chosen  to 

point  in  the  X  direction.  Therefore,  the  general  formulation  of  the  line  of  the  node  is 
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N  = 


cos(Q) 

sin(Q) 

0 


1 

0  if  z  =  0  a  y  =  0 

0 

p  sin  j  cos  %  +  z  (cos  t/cos^sin(s7/L)  +  sin4rcos(s7/L)) 
p  sin  jsin<^  +  z(cos  jsin^rsin(s7/L)-cos^cos(s7/L))  otherwise 
0 


(3.9) 


Note  that,  if  the  NKO  orbital  plane  is  parallel  to  the  equator  and  z  ^  0 ,  Eq.  (3.9)  becomes  the  following. 


N(j  =  0)=, 


cos  ^  sin  (mt)  +  sin  E,  cos  (mt ) 
-cos  E,  cos  (mt)  +  sin  E,  sin  ( mt ) 
0 


(3.10) 


Because  the  NKO  RAAN  E,  is  not  defined  for  j  =  0 ,  it  can  be  arbitrarily  set  to  zero  and  Eq.  (3.10)  becomes  the 
following,  which  is  the  description  already  obtained  in  the  simplified  case  (Eq.  (2.10)). 


n„  =  A 


sin  (mt) 
-cos  (mt) 
0 


(3.11) 


The  right  ascension  of  the  ascending  node  is  simply  given  by 

Q  =  arctan  (n-X,N-y)  = 

in  which 


0  if  z  =  ()  a]  =  0 
Cl  otherwise 


(3.12) 


Cl*  =arctan(psin  ycos^  +  ^cos  j  cos  £  sin  (art)  +  sin  £  cos  (nrt)),p  sin  jsin£  +  z(cos  j  sin  ^  sin  (mt)- cos  ^  cos  (mt))j 

(3.13) 

Note  that  the  term  1/N  has  been  neglected  in  Eq.  (3.13)  because  it  is  always  a  positive  value. 

Note  that,  if  z  =  0  but  j  ^  0 ,  the  osculating  RAAN  equals  that  of  the  NKO: 


Cl(z  =  0)  =  £ 


(3.14) 


On  the  other  hand,  if  the  NKO  orbital  plane  is  parallel  to  the  equator,  Eq.  (3.12)  becomes  Eq.  (3.15),  which  is  the 
description  already  obtained  in  the  simplified  case  (Eq.  (2.23)). 


n(j  =  £  =  o) 


0  ifz  =  0 

arctan  ( sin  (mt)  sgn  ( z ) ,  -  cos  (mt)  sgn  ( z ))  otherwise 


(3.15) 


Considering  that  it  is  arbitrarily  chosen  that  E,  =  0  for  a  NKO  with  the  orbital  plane  parallel  to  the  equatorial  plane, 
Eq.  (3.12)  could  simply  be  rewritten  as  Cl  =  Cl  ,  taking  into  account  Eq.  (3.14).  However,  it  can  be  seen  that 

Q*  (j  =  z  =  0)  =  arctan  (0,0)  (3.16) 

Numerically,  Eq.  (3.16)  often  assumes  the  zero  value  but  the  expression  in  Eq.  (3.16)  is  indefinite,  if  solved  analyti¬ 
cally.  Moreover,  if  one  wants  to  compute  the  limits  to  find  the  exact  analytic  value  of  Eq.  (3.16),  it  can  be  demonstrated 
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that  the  limits  are  different  depending  on  the  order  of  the  limits  themselves.  If  the  previous  statement  is  true,  the 
following  holds. 

limflim^Y  W  limflimQ*)  (3.17) 

j— >0  \  z— >0  /  z— >0  \  j—>0  ] 

If  Eq.  (3.17)  is  true,  it  means  that  Q*  ( j  =  z  =  0)  =  arctan(0,0)  is  indeterminate  and,  therefore,  the  general  analyt¬ 
ical  expression  for  the  osculating  right  ascension  of  the  ascending  node  is  the  one  given  in  Eq.  (3.12). 

The  left-hand  side  of  Eq.  (3.17)  is 

limllimQ  )  =  limfarctanfpsin jcos^psin jsin^))^^  (3.18) 

Note  that  the  second  limit  in  Eq.  (3.18)  has  been  solved  using  the  L’Hopital’s  rule.  The  right-hand  side  of  Eq.  (3.17) 
is 

lim  |limQ*  j  =  lim  ^arctan^(cos^sin(^)  +  sin<^cos(^)),  z(W^sin(s7f)-cos^cos(W)))j  = 

=  arctan(cos^sin(s7r)  +  sin4rcos(cjf),sin4:sin(s7r)-cos^cos(s7r))  =  (3.19) 

=  arctan(sin(ft7t  +  £),-cos(©?  +  £))  =  arctan(cos(Wf  +  %  - 7t / 2) ,sin(mt  +  £-;r/2))  => 


lim ( limQ  \  =  mt  +  <E-— ^<E  =  lim(  limQ*  ) 

z— »0+  \  j-> 0  /  2  J-^0  V  <  J 


(3.20) 


Therefore,  Eq.  (3.20)  demonstrates  that  the  general  analytical  expression  for  the  osculating  right  ascension  of  the 
ascending  node  is  the  one  given  in  Eq.  (3.12). 

Lastly,  the  term  Q*  can  be  rewritten  as 


Jarccos  (NjN)  ifNy>0 

\2ir  -  arccos  ( Nx  /  N )  otherwise 


(3.21) 


in  which 

N  =  <Jz2  cos2  (s7f)  +  (/?sin  j  +  zcos  jsin(sj^))2 
<  Nx  =  p  sin  jcos^  +  z^cos  j  cos  d;  sin  (eat)  +  sin  £  cos  (orf))  (3.22) 

Ny  =  p  sin  jsin^  +  z(cos  y  sin^sin(G7t)-cos<f:cos(G7t)) 


3.  Eccentricity  ( e ) 

The  eccentricity  vector  is  given  by 


A  vxh 

e  =  -r- 1 - 

M 


(3.23) 


Using  the  vector  triple  product  rule  and  remembering  that  vr  =  0 ,  the  expression  of  the  eccentricity  vector  in  Eq. 
(3.23)  can  be  rewritten  as 
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V 


_ \ 

M  r 


_2  2 

m  p 


\ 

1 

p  (cos  <!;  cos  (mt)- cos  jsin4rsin(c7^))  + zsin  y'sin^ 

>o(sin4:cos(^)  +  cos  y'cos^sinfW^-zsin  j  cos^ 

sIP  +z  ) 

p  sin  y'sin(67t)  +  zcos  j 

(3.24) 


Therefore,  the  value  of  the  eccentricity  is  given  by  the  norm  of  the  eccentricity  vector,  as 


e  = 


p-m2p2jprTz2 

P 


(3.25) 


Note  that  a  circular  planar  Keplerian  orbit  is  characterized  by  p  =  zu2 p3 ,  which  satisfies  the  zero-eccentricity 
condition.  However,  the  eccentricity  vector  is  not  defined  in  this  case  and  this  is  an  issue  in  the  following  computation 
of  the  true  anomaly.  Therefore,  the  eccentricity  unit  vector,  in  this  case,  is  arbitrarily  chosen  to  point  in  the  X  direc¬ 
tion.  Therefore,  the  general  formulation  of  the  eccentricity  vector  is 


e 


1 

0 

0 


2 

zu  p 


\ 

1 


/?(cos<^cos(67t)-cos  jsin^sin(s7/L))  +  zsin  j  sin^ 
/?(sin<^cos(67t)  +  cos  y  cos<^sin(67t))-zsin  j  cos^ 
p  sin  y  sin(W)  +  zcos  j 


if  z  =  0  a//  =  m2p3 


otherwise 


(3.26) 


4.  Semi-major  axis  (a) 

The  semi -major  axis  is  given  by 


H2 

m\Ip2+z2 

//(i-e2) 

2  ju-m2p2s]p2  +z2 

(3.27) 


5.  Argument  of  perigee  ( co ) 

The  argument  of  perigee  is  given  by 


arccos^A  -ej 
2 n - arccos \N  -e 


if  ez  >  0 
otherwise 


(3.28) 


After  algebraic  manipulations  and  arbitrarily  choosing  co  =  0  in  the  case  of  osculating  circular  orbit,  the  general 
expression  for  the  argument  of  perigee  is 


co  = 


0 

CO 

In -co* 


if  p  =  TU2  p2  ffFTz2 

if  [^p-w2  p2  fp2  +  z2  j(zcos  7  +  /7sin  j'sin(sjf))<0 
otherwise 


(3.29) 


in  which 
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co  =  arccos 


sin  j  cos  (zut)^p2  +z2 


^z2  cos2  (zat) +  (p sin  j  +  zcos  jsin(zat)J 


rsgn 


(3.30) 


Note  that,  if  the  NKO  orbital  plane  is  parallel  to  the  equator  and  z  ^  0 ,  Eq.  (3.29)  becomes  the  following,  which 
is  the  description  already  obtained  in  the  simplified  case  (Eq.  (2.23)). 


CD 


0= °)= 


nj  2  if  >  0  a//  <  m2p2yjp2  +  z2  <  0a  p>  m2p2yjp2  +z2  J 

3tt/2  if  >  0  a//  >  m2p2yjp2  +  z2  j  <  0  a//  <  m2p2yjp2  +z2  J 


(3.31) 


If  the  NKO  orbital  plane  is  parallel  to  the  equator  and  z  =  0,  Eq.  (3.29)  becomes  the  following,  which  is  the 
description  already  obtained  in  the  simplified  case  (Eq.  (2.23)). 


liml  lim co)  ■ 

y— >o  \  z— »o  / 


0  if  p  -  m2 p3 

cut  if  p  <  zu2 p3 

71  +  tnt  if  p>zu  2p3 


(3.32) 


Note  that  Eq.  (3.31)  has  been  obtained  by  simply  substituting  j  =  0  within  Eq.  (3.29).  On  the  other  hand,  Eq. 

(3.32)  has  been  obtained  by  computing  two  limits.  Therefore,  the  general  complete  expression  for  the  argument  of 
perigee  is  given  by  unifying  Eqs.  (3.29)  and  (3.32),  as  follows. 


co  = 


0 

if  p  =  zu2  p 2^p2  +  z2 

cut 

ifj  =  0Az  =  0/\p<  zu2 p3 

71  +  VJt 

lfj  =  0AZ=0Ap>  zu2  p3 

CD* 

if  +Z2^[z ( 

17T-CD* 

otherwise 

(3.33) 


6.  True  anomaly  (  3 ) 

The  true  anomaly  is  given  by 


cos 


(*)= 


e  r 


(3.34) 


In  the  case  of  a  circular  planar  Keplerian  orbit,  Eq.  (3.34)  becomes  the  following.  This  is  due  to  the  choice  of  the 
eccentricity  vector  shown  in  Eq.  (3.26). 


cos(l9)  =  cos  {put )  =>  3  =  cut 

On  the  other  hand,  a  NKO  is  characterized  by  the  expression  of  the  true  anomaly  given  as 

3  =  arccos  sgn  {^p  -  zu2  p2  y] p 2  +z2  jj 

Eq.  (3.36)  becomes 


3  = 


\  0  if  p<m2p2yjp2  +z2 
[tT  if  p  >  ZU2 p2  y ]p 2  +  z 


(3.35) 


(3.36) 


(3.37) 
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Note  that  the  case  |z^0a//  =  m2p2«Jp2  +z2  j  is  not  taken  into  account  inside  Eq.  (3.37).  From  Eq.  (3.36),  it  can 

be  seen  that  i9|z^0a//  =  m2  p2  *Jp2  +z2^j  =  7tj  2 .  However,  it  does  not  exist  any  case  for  which  p  -  zu2  p2  ^p2  +  z 

if  z  ^  0 ,  i.e.  there  are  no  osculating  circular  orbits  for  non-Keplerian  orbits.  Therefore,  the  general  expression  for  the 
true  anomaly  can  be  rewritten  by  unifying  Eqs.  (3.35)  and  (3.37)  as 


mt 

0 

71 


if  z  =  0  a//  =  m2p3 
if  p<m2  p2  ^p2  +z2 
if  p  >  m2p2^p2  +  z2 


(3.38) 
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Summary 

The  generalized  forward  map  from  NKOs  to  osculating  classical  Keplerian  elements  described  in  Section  3.1.1  is 
summarized  as: 


a  = 


+  z2 


2p-m2  p2  \/p2  +z: 
p-zu2  p2*(pAz2 


i  -  arccos 


p cos  j  -  z sin  ysin(W) 


p2+z2 


Q: 


CO  - 


*9  = 


0  if  z  =  0  a  j  =  0 
Q*  otherwise 

if  p  =  m2p 2Pp2  +  z2 

if/  =  0Az  =  0A//<  m2  p3 
ifj-0Az  =  0Ap>  m2p3 

if  ^ju - m2 p2 <sjp2  +  z2  j(zcos  j  +  p sin  7'sin(^))<0 
otherwise 
mt  if  z  =  0  a  p  =  m2  p3 
0  if  jU  <  m2p2^Jp2  +  z2 


0 

mt 

7r  +  mt 
GO 

In- co 


n 


if  p  >  m2 p2 ^[pAz2 


(3.39) 


in  which 


'  „  farccos(iVx/iV)  ifA^>0 
|2;r  -  arccos  /Af )  otherwise 

Af  =  ^ jz 2  cos2  (ft7t)  +  (/?sin  7  +  zcos  7'sin(£7t))2 
A^  =  /?sin  j  cos  <^  +  z  (cos  jcos^sin^/^  +  sin^cosjW)) 
'  Ny  =  /9sin  j  sin  <^  +  z  (cos  j  sin  £  sin  {i&t)-  cos  £  cos  (cjt)) 


00  -  arccos 


7  = 


(-r  sgn  [p  -  ™ 2  p2 

sin  j  cos  (mt)^ 


/P2  +  Z2 


ijz2  cos2  (cj^  +  ^sin  7 +  zcos7sin(cj/L))2 


(3.40) 
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3.1.2.  Modified  Equinoctial  Elements 

Starting  from  the  expressions  of  the  Cartesian  position  and  velocity  vectors  (Eqs.  (3.2)  -  (3.3)),  the  orbital  angular 
momentum  (Eq.  (3.4))  and  the  analytical  map  of  the  osculating  classical  Keplerian  elements  described  in  Section  3.1.1 
(Eq.  (3.39)),  a  map  of  the  osculating  MEE  from  Cartesian  position  and  velocity  is  derived. 

1.  Semi-latus  rectum  (  p ) 


|ftf  _»VV+£) 

M  M 


(3.41) 


2.  In-plane  element  (  f  ) 


In  order  to  derive  the  general  formulation  for  the  in-plane  element  /  ,  let  us  first  consider  the  case  of  zero-dis¬ 
placement  orbits  (i.e.  z  =  0 ). 

From  the  formulation  of  Q  shown  in  Eq.  (3.39),  the  value  of  the  right  ascension  of  the  ascending  node  for  a  zero- 
displacement  orbit  is  simply 


(3.42) 


From  the  formulation  of  co  shown  in  Eq.  (3.39)  and  noting  that  p sin  j  >  0  always,  the  value  of  the  argument  of 
perigee  for  a  zero-displacement  orbit  is 


= 


td t  if  ju  <  m2  p3 


\jr  +  tnt  if  ju>tn2p3 
From  Eqs.  (3.42)  -  (3.43),  the  value  of  the  longitude  of  perigee  [ta  +  Q]  is 


(3.43) 


[®  +  nL=o  - 


%  +  UJt  if// <67 2  p3 

%  +  7T  +  znt  if  //  >  zu2  p3 


(3.44) 


cos(&>  +  Q)|  _q  =  < 


COs(£  +  G7f)  if  p  <  VJ2 p3 
-cos(<^  +  67f)  if  ju>zn2p3 


(3.45) 


From  Eq.  (3.45)  and  remembering  the  expression  of  the  eccentricity  shown  in  Eq.  (3.39),  the  expression  of  the  in¬ 
plane  element  /  for  a  zero-displacement  orbit  is 


f\z= 0  =[eCOS(®  +  Q)l=0  =~ - -C0S(/  +  B7 1)  : 

P 

2  3  _ 

/ 1  _0  =  m  P — —  [cos  ^cos^/)  -sin  sin  (/at)] 
// 


(3.46) 


(3.47) 


If  the  NKO  orbital  plane  is  parallel  to  the  equator  (i.e.  j  =  0  but  also  £,=()),  the  same  procedure  described  in 
Section  2. 1.1. 2  can  be  used  to  obtain 


/L 


uj2 p2^p2  +z2  -// 


cos 


(m) 


(3.48) 
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Note  that,  if  the  orbital  plane  is  parallel  to  the  equator  and  z  =  0 ,  Eqs.  (3.47)  and  (3.48)  are  equivalent. 

Let  us  consider  now  a  displaced  orbit.  From  the  formulation  of  Q  shown  in  Eq.  (3.39),  the  value  of  the  right 
ascension  of  the  ascending  node  for  a  displaced  orbit  is 

QL=Q*  (3-49> 

From  the  formulation  of  co  shown  in  Eq.  (3.39),  the  value  of  the  argument  of  perigee  for  a  displaced  orbit  is 


co\ 


z*  0 


arccos(^)  if  p-m2p2^p2  +z  <0azcos  j  +  p sin  jsm[rut)  >  0 

arccos  (-y)  if  p  -  m2  p2  «Jp2  +z2  >  0  a  z  cos  j  +  p  sin  j  sin  (rut)  <  0 

-  arccos  (-y)  if  p  -  vj2  p2  ^p2  +z2  >0az  cos  j  +  p  sin  j  sin  [rut )  >  0 

-  arccos  (x)  if  p-vj2p2^p2  +z  <0azcos  j  +  p  sin  jsm(mt)  <  0 


(3.50) 


in  which  y  has  been  defined  in  Eq.  (3.40). 

From  Eqs.  (3.49)  -  (3.50)  and  recalling  that  arccos  (— jc)  =  ;r- arccos  (v)  ,  the  value  of  the  cosine  of  the  longitude 
of  perigee  cos(&>  +  Q)|  is 


cos(&>  +  Q)|  ^ 


cos  +  arccos  (/))  if  p  -  m2 p2  ^jp2  +  z2  <  0  a  z  cos  j  +  p  sin  ysin(W)>0 
-cos  (Q*  -arccos  (/))  if  p  —  m2p2  ^p2  +  z2  >  0  a  z  cos  j  +  p  sin  ysin(W)<0 
-  cos  (Q*  +  arccos  (/))  if  p  -  m2p2  ^p2  +  z2  >  0  a  z  cos  j  +  p  sin  ysin(W)>0 
cos  (Q*  -  arccos  (/))  if  p  -  m2 p2  4+  +  z2  <  0  a  z  cos  j  +  p  sin  j  sin  (jut)  <  0 


(3.51) 


From  Eq.  (3.51)  and  remembering  the  expression  of  the  eccentricity  shown  in  Eq.  (3.39),  the  expression  of  the  in¬ 
plane  element  /  for  a  displaced  orbit  is 


cos  (Q*  +  arccos  (r))  if  z  cos  j  +  p  sin  j sin( rut)  >0 

;  ,  (3.52) 

cos  ( Q*  -arccos  (7)  j  if  z  cos  j  +  p  sin  jsin(rut)<0 

Note  that,  if  z  =  0 ,  Eq.  (3.52)  becomes  Eq.  (3.47).  On  the  other  hand,  if  j  =  0 ,  Eq.  (3.52)  becomes  Eq.  (3.48),  which 
is  the  expression  found  in  the  simplified  case  (Eq.  (2.76)).  This  can  be  easily  demonstrated  as  follows. 

Let  us  first  consider  the  z  +  0  case.  It  can  be  seen  that  y |  =  0  .  That  is,  Eq.  (3.52)  becomes 


fLo=- 


m"p2_jpf+zf_  -ju 

P 


4- 


<n2p24p2  +  z2  -p 


COS 

[n-+4| 

1 

l  2) 

1 

f  t 

cos 

1 

l  2  J 

if  z  >0 
if  z  <  0 


(3.53) 


As  demonstrated  in  Eq.  (3.19),  Q*|  q  =  %  +  mt  —  nl2  if  z  >  0 .  Following  the  same  procedure,  it  can  be  demon¬ 
strated  that  the  following  holds. 


Q* 


\j=0 


(<^  +  mt-7r/2  ifz>0 
y^  +  vut  +  n/ 2  if  z  <  0 


(3.54) 
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Substituting  Eq.  (3.54)  into  Eq.  (3.53)  and  recalling  that  it  is  arbitrarily  chosen  that  E,  =  0  if  j  =  0 ,  the  expression 
of  the  in-plane  element  /  becomes  the  same  as  that  shown  in  Eq.  (3.48). 

Let  us  consider  now  the  z  =  0  case.  It  can  be  seen  that  y\  Q  =  cos (gtt) .  Recalling  that  Q|  _^_o  =  0  (Eq.  (3.39)), 
Eq.  (3.52)  becomes  Eq.  (3.48). 

Therefore,  the  general  expression  for  the  in-plane  element  /  is 

cos  {tut)  if  z  =  0  a  j  -  0 

cos^Q  +  arccos  (;k))  if  zcos  j  +  psm  >  0  (3.55) 

cos^Q*  -  arccos  (;k))  if  z  cos  j  +  p sin  jsm(mt)<0 

It  is  important  to  note  that  Eq.  (3.55)  is  a  piecewise  formulation  characterized  by  nested  direct  and  inverse  trigo¬ 
nometric  functions.  Moreover,  Q*  itself  is  defined  through  a  piecewise  law  (Eq.  (3.40)).  For  these  reasons,  an  alter¬ 
native  formulation  can  be  found,  after  some  mathematical  manipulations,  as  follows. 

*  y 

Using  the  same  nomenclature  used  for  Q  ,  /  can  be  rewritten  as  /  =  —  .  That  is, 

N 

y-  sin  j  cos  (zut)*Jp2  +z2  (3.56) 

Taking  into  account  the  second  condition  in  Eq.  (3.55)  and  the  condition  Ny  >  0 ,  the  expression  for  the  in-plane 
element  /  is 


f=- 


m^p2^p^  +  z"  -  p 

P 


f  = 


nr2  p2  *Jp2  +  z  -p 
P 


-cos 


/ 

(  Nx) 

f  y\] 

arccos 

— 

+  arccos  — 

V 

{N  ) 

UJJ 

(3.57) 


Let  us  recall  some  useful  relations  between  trigonometric  functions. 


arccos  (x)  =  2  tan  1 


tan  1(x)±tan  1(y)  =  tan  1 


2 

cos  ^2  tan1  (x))  =  — — - 

Using  Eqs.  (3.58)-(3.59)  and  after  some  algebraic  manipulations,  the  following  relation  can  be  obtained. 

f  „  „  A 


Ji  —  x 

(3.58) 

1  +  X 

) 

rlf  x±u 

U+^yJ 

(3.59) 

2,-. 
+  1 

(3.60) 

arccos 


N 


+  arccos 


f-j^l  =  ...  =  2  tan-1 


lN2-N2(N  +  y)  +  yl 

’In2-/2 

(N  +  Nx) 

(N  +  Nx)(N  +  r)-J 

\n2-nx2) 

(n2-/2) 

(3.61) 


Applying  Eq.  (3.60)  to  Eq.  (3.61),  the  following  can  be  obtained. 

Nxr-J{N2~Nx2)(N2-?) 


cos 

( 

arccos 

fNx 

)  (  y  ^ 

V 

)  UJJ 

N1 


(3.62) 
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Following  the  same  procedure  described  in  Eqs.  (3.57)  -  (3.62),  the  formulation  of  the  in-plane  element  /  can  be 
rewritten  as 


f  _m2p14p2  +  Z2  -p 
1  pN2 


N 2  cos 

(m) 

if  z  —  0  a  j  -  0 

J{N2~N*2) 

(N2-f) 

if  -(zcos  7  +  /?  sin  j  sin(mt)^)>0 

(3.63) 

nj  +  . 

(N2-f2) 

if  -(zcos  7  +  /?  sin  j sin(mt)^)  <0 

The  formulation  of  the  in-plane  MEE  /  shown  in  Eq.  (3.63)  is  free  from  any  nested  direct  and  inverse  trigono¬ 
metric  function.  Nevertheless,  this  formulation  is  not  any  simpler  than  the  one  shown  in  Eq.  (3.55).  Moreover,  from  a 
numerical  point  of  view,  the  computation  of  the  latter  is  faster.  For  these  reasons,  Eq.  (3.55)  is  considered  as  the  final 
formulation  of  the  in-plane  MEE  /  . 


3.  In-plane  element  (  g  ) 

As  for  the  case  of  the  in-plane  element  /  ,  let  us  derive  the  general  mapping  for  the  in-plane  element  g  consid¬ 
ering  the  in-plane  and  out-of-plane  orbits  separately  at  first.  In  order  to  derive  the  general  formulation  for  the  in-plane 
element  g  ,  let  us  first  consider  the  case  of  zero-displacement  orbits  (i.e.  z  =  0 ). 


From  Eq.  (3.44),  the  value  of  the  longitude  of  perigee  [&>  +  Q]  Q  is 

r  .  f  E  +  mt  if  u<zu2p3 

Z  [Z  +  TT  +  m  if  JU>ZU2p 


/  v,  I  sin (E  +  mt)  if  u<xn2p 

sin(^  +  Q)  =  H  H 

v  y|‘=°  |-sin(^  +  ^)  ifp 


2  3 

>  m  p 


(3.64) 


(3.65) 


From  Eq.  (3.65)  and  remembering  the  expression  of  the  eccentricity  shown  in  Eq.  (3.39),  the  expression  of  the  in¬ 
plane  element  g  for  a  zero-displacement  orbit  is 


U=[esin(®+Q)Lo=  P/u  //sin(^  +  S7f)  => 

(3.66) 

2  3 

#L_0  =  — ^ — —  [sin£cos(E7t)  +  cos£sin(G7t)] 

P 

(3.67) 

If  the  NKO  orbital  plane  is  parallel  to  the  equator  (i.e.  j  =  0  but  also  £  =  0 ),  the  same  procedure  described  in 
Section  2. 1.1. 2  can  be  used  to  obtain 


,  m 2p1Ip2Tz2~m  .  .  s 

gU  -  , - -sm M 


(3.68) 


Note  that,  if  the  orbital  plane  is  parallel  to  the  equator  and  z  =  0 ,  Eqs.  (3.67)  and  (3.68)  are  equivalent. 

Let  us  consider  now  a  displaced  orbit.  Following  what  have  been  done  for  the  case  of  the  in -plane  element  /  ,  the 
value  of  the  sine  of  the  longitude  of  perigee  sin  (  &>  +  Q)|  is 
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sin(®  +  Q)L, 


sin  (Q*  +  arccos  (7))  if  p-  + p 2jp2  +  Z2  <0azcos  j  +  p  sin  j  sin(mt)>0 
-sin(Q*  -  arccos  (;r))  if//-  m 2p2fp2++  >0  a  z  cos  j  +  p  sin  j  sin(W)  <0 
-sin(Q*  +  arccos  (;r))  if  p-  m2p2jp+72  >0azcos  j  +  p  sin  j  sin(W)>0 
sin^Q*  -  arccos  (x))  if  p-eu1  p2  \[p2  +  z2  <0azcos  j  +  psin  jsm(mt)<0 


(3.69) 


From  Eq.  (3.69)  and  remembering  the  expression  of  the  eccentricity  shown  in  Eq.  (3.39),  the  expression  of  the  in¬ 
plane  element  g  for  a  displaced  orbit  is 


sin(Q*  +  arccos  (/))  if  z  cos  j  +  p  sin  jsm(mt)>0 

*  (3-7°) 
sin  ( Q* -arccos  (f)J  if  z  cos  j  +  p  sin  y'sin(£7t)<0 

Note  that,  if  z  =  0 ,  Eq.  (3.70)  becomes  Eq.  (3.67).  On  the  other  hand,  if  j  =  0 ,  Eq.  (3.70)  becomes  Eq.  (3.68),  which 
was  found  in  the  simplified  case  (Eq.  (2.76)).  The  demonstration  of  this  follows  the  same  steps  shown  in  the  case  of 
the  in-plane  element  /  . 

Therefore,  the  general  expression  for  the  in-plane  element  g  is 


*L=- 


^p2jp^rZ_  ~p  _ 

M 


sin  {tut)  if  z  =  0  a  j  =  0 

sin^Q*  +  arccos  (f))  if  zcosj  +  psinjsm(mt)  >  0  (3.71) 

sin  (Q*  -  arccos  (/))  if  z  cos  j  +  p  sin  j  sin  (mt )  <  0 

As  for  the  in-plane  element  /  ,  it  can  be  noted  that  Eq.  (3.55)  is  a  piecewise  formulation  characterized  by  nested 
direct  and  inverse  trigonometric  functions.  Nevertheless,  it  can  be  demonstrated  that  a  procedure  similar  to  that  for 
the  in-plane  element  /  does  not  simplify  the  formulation.  Therefore,  such  a  mathematical  derivation  is  not  shown 
here  for  the  sake  of  conciseness  and  Eq.  (3.71)  is  considered  as  the  final  formulation  for  the  in-plane  MEE  g  . 


zu^p2^p^  +  z^  -  P 


4.  Out-of-plane  element  (  h ) 

Before  deriving  the  expression  for  the  out-of-plane  element  h  ,  let  us  recall  the  following  identity. 


tan 


l-cos(x) 
V2  J  sin(x) 


Recalling  the  formulation  of  the  inclination  shown  in  Eq.  (3.39)  as  i  =  arccos 
pression  of  tan(//2)  is,  therefore,  given  by 


p  cos  j-zsin  jsin(mt) 


■j, 


p1  +  z 


(3.72) 


,  the  ex- 
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p  cos  y'-zsin  jsin(nrt) 


“(i)- 


j?+i 


r  .  .  . .  / 

p  cos  y-zsin  j  sm  ymt) 


sin 


arccos 


4p2 


(3.73) 


JJ 


^  p  cos  y'-zsin  ysin(sj/L) 


\[p2  +z2 


P2  +  Z2 


p2  +  z2  ~(p cos  j  —  z  sin  ysin(c7t)) 


tan 


( i\  4p2  +z2  ~{pcos  j~z  sin 

-^p2  +  z2  ~(p cos y'-zsin ysin(c7?))2 


(3.74) 


Noting  that  the  denominator  of  Eq.  (3.74)  can  be  rewritten  as  4 A2  -B2  =  ^(A  +  Z?)(A-Z?)  ,  Eq.  (3.74)  can  be 
simplified  as 


i  |  |^/y02+z2  (/^  cos  y'-zsin  ysin(cn)) 


Up2  +  z2  +(>ocos  y'-zsin  ysin(^)) 
From  Eq.  (3.75)  and  recalling  the  value  of 


(3.75) 


cos 


(Q): 


1  if  z  =  0  a  y  =  0 

/?sin  ycos^  +  z(cos  y'cos^sin^^  +  sin^cos^r)) 


(3.76) 


sjz2  cos2  (W )  +  (/?  sin  y  +  zcos  y  sin (W))2 
introduced  in  Eq.  (3.9),  the  expression  of  the  out-of-plane  element  /z  is 


otherwise 


h  = 


~(p 


cos  y  -  z  sin  y  sin 


M) 


'x//?2  +  z2  + 


cos  y-zsm 


./'sin  (art)) 


if  z  =  0  a  y  =  0 


yosin  ycos£  +  z(cos  y'cos^sin(^^  +  sin^cos(67t)) 


^ jz 2  cos2  (W)  +  (/?sin  y'  +  zcos  y'sin(sjr))2 


otherwise 


(3.77) 


Note  that  the  piecewise  formulation  in  Eq.  (3.77)  can  be  removed  by  performing  the  two-variable  limit,  as  follows. 
Let  us  consider  only  the  second  branch  of  Eq.  (3.77)  to  demonstrate  that 


lim  h  =  0  (3.78) 

(*.j)-K°.  o) 


First,  let  us  rewrite  the  two  variables  z,  j  in  polar  coordinates  as 


jz  =  r cos  3 
[j  =  r  sin  8 

Using  Eq.  (3.79),  the  second  branch  of  Eq.  (3.77)  becomes 


(3.79) 
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h  = 


E 

jp2  +  r 2  cos  *9  -  p  cos  | 

(r  sin^j 

+  rcosi9sin| 

(r  sin^j 

sin(^) 

p2  +  r2  cos2  *9  +  /?cos 

^rsin»9 

)-rcost9sin 

sin»9 

)sin(67t) 

yocos^sinj 

(rsin#) 

|  +  rcos^cos^cos 

(r  sin^j 

|  sin  (mt)  +  sin  £  cos  [zut )  j 

i 

lr2  cos2  3 cos2  (m)  +  i 

^sin^f  sin*9j 

+  f  cosi9cos| 

(r  sin^j 

|sin(67t)j 

The  limit  shown  in  Eq.  (3.78)  can  be  now  re-written,  using  the  transformations  in  Eq.  (3.79),  as 

lim  h  =  0 

r— »0+ 


(3.80) 


(3.81) 


It  can  be  easily  demonstrated  that  Eq.  (3.81)  is  verified  for  any  value  of  3  and,  therefore,  the  two-variable  limit 
shown  in  Eq.  (3.78)  is  verified  as  well. 

Considering  that  ^  hin  ^  h  -  0 ,  the  expression  of  the  out-of-plane  element  h  is  simply 


k 

fp2+z2  —  (p cos 7-zsin 7'sin(W))  p sin ^cos^-l- z(cos 7COS^sin(£7/^  +  sin^cos(£7t)) 

H 

^ p2  +  z  +(yOCOS7 -zsin  jsinipnt^  ^ 

1 z2  cos2  (c7^)  +  (yC>sin  7  +  zcos  7'sin(G7t))2 

(3.82) 


5.  Out-of-plane  element  (  k  ) 

As  computed  in  Eq.  (3.75)  for  the  out-of-plane  element  h  ,  the  expression  of  tan(//2)  is 


tan 


np2  +^2  cos  j  —  z  sin  j  sin  )) 

v  ^jp2  +  z2  +(p cos  7-zsin  jsin(^)) 


From  Eq.  (3.83)  and  recalling  the  value  of 


sin(Q)  = 


0  if  z  =  0  a  7  =  0 

p  sin  7'sin4r  +  z(cos  j  sin  ^  sin  (gt?)-  cos  %  cos  (W)) 


<Jz2  cos2  (mt)  +  (p sin  j  +  zcos  j  sin  (mt^ 
introduced  in  Eq.  (3.9),  the  expression  of  the  out-of-plane  element  k  is 


otherwise 


(3.83) 


(3.84) 


k  = 


l 


p2  +  z2  ~(p  cos  j  —  z  sin  j  sin 


M) 


y  \[p2  +  Z2  +(pcos  j-z sin 7 sin (arf  )) 


if  z  =  0  a  y'  =  0 


/?sin  7‘sin4r  +  z(cos  j  sin  sin  (tat)-  cos  cos  (art)) 


v 


otherwise 


z2  cos2  (^)  +  (/?sin  7  +  z  cos  7'sin(sjr)) 


(3.85) 
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Note  that  the  piecewise  formulation  in  Eq.  (3.85)  can  be  removed  by  performing  the  two-variable  limit,  as  demon¬ 
strated  for  the  out-of-plane  element  h  .  Let  us  consider  only  the  second  branch  of  Eq.  (3.85)  to  demonstrate  that 


lim  k  =  0 

M-K  °’0) 

First,  let  us  rewrite  the  two  variables  z,j  in  polar  coordinates  as 

z  =  r cos  3 
j  =  r  sin  3 

Using  Eq.  (3.87),  the  second  branch  of  Eq.  (3.85)  becomes 


k  - 


A 

jp2  +  r 2  cos  3  -  p  cos  ( 

rsin^j 

+  rcosz9sin| 

[r  sin^j 

sin  (mt) 

p2  +  r2  cos2  3  +  p  cos 

[r  sini9 

)-rcosz9sin 

(r  sini9 

)sin(W) 

p  sin 


E,  sin  (r  sin  3  j  +  f  cos  3  ^sin  E,  cos  {r  sin  3  j  sin  ( mt )  -  cos  E,  cos  (mt )  j 


Jr2  cos2  i9cos2  (art) +  (p  sin  (r  sin^j-t  r  cos^cos^r  sinl9jsin(s7/L)j 
The  limit  shown  in  Eq.  (3.86)  can  now  be  re-written,  using  the  transformations  in  Eq.  (3.87),  as 


lim  k  =  0 

r— »0+ 


(3.86) 


(3.87) 


(3.88) 


(3.89) 


It  can  be  easily  demonstrated  that  Eq.  (3.89)  is  verified  for  any  value  of  3  and,  therefore,  the  two-variable  limit 
shown  in  Eq.  (3.86)  is  verified  as  well. 

Considering  that  ^  lim  ^  k  =  0  ,  the  expression  of  the  out-of-plane  element  k  is  simply 


k  Up2  +  z2  —(p  cos  j-zsin  ysin(zzTt))  psin  j  sin  £  +  z(  cos  j  sin  %sin(znt )- cos  E,  cos  (cTt))  ^399^ 

I  *Jp2  +z2  +(pcosj-zsinjsin(m))  L2  cos2  (mt)  +  (p sin  j  +  z cos  y  sin  (art))2 


6.  True  longitude  ( L ) 

In  order  to  derive  the  general  formulation  for  the  true  longitude  L  ,  let  us  consider  the  two  cases  z  =  0  and  z  ^  0 
,  respectively. 

From  the  formulations  of  Q  and  00  shown  in  Eq.  (3.39),  the  value  of  the  longitude  of  perigee  [&>  +  Q]  Q  is 


[®+QL 


E,  Up-  w2 p3 
<  %  +  znt  if  p  <  m2 p3 

%  +  7t  +  mt  iip>zu2pi 


(3.91) 


Therefore,  the  value  of  the  true  longitude  in  the  case  of  zero-displacement  NKO  is  obtained  from  Eqs.  (3.39)  and 
(3.91)  as 

L\z=Q=£  +  m  (3.92) 
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Let  us  consider  now  a  displaced  orbit.  From  Eqs.  (3.49),  (3.50),  and  (3.39)  and  recalling  that 
arccos  (-v)  =  tt- arccos  (v) ,  the  value  of  the  true  longitude  for  an  out-of-plane  displaced  orbit  is 


4,»  = 


j  Q  +  arccos  (y)  if  z  cos  j  +  p  sin  j  sin  (mt )  >  0 
[Q*  -arccos^)  if  z cos  j  +  p sin  jsin(mt)<0 


(3.93) 


It  is  possible  to  demonstrate  that  Eq.  (3.93)  becomes  Eq.  (3.92)  in  the  case  of  zero-displacement  orbit.  However, 
the  case  of  zero-displacement  NKOs  with  the  orbital  plane  parallel  to  the  equator  must  be  treated  in  a  different  way, 
as  demonstrated  in  the  calculation  of  the  right  ascension  of  the  ascending  node  and  shown  in  Eq.  (3.39).  Therefore, 
the  general  expression  for  the  true  longitude  is 


L  =  \ 


mt  if  z  =  0  a  j  -  0 

Q  +  arccos  (^)  if  z  cos  j  +  p  sin  j  sin  ( mt )  >  0 
Q*  -  arccos  (y)  if  z  cos  j  +  p  sin  j  sin  (mt)  <  0 


(3.94) 


An  alternative  formulation  for  the  true  longitude  can  be  obtained  which  does  not  involve  inverse  trigonometric 
functions,  as  described  in  the  case  of  the  in-plane  element  /  .  However,  as  discussed  for  the  in-plane  modified  equi¬ 
noctial  elements,  it  can  be  demonstrated  that  such  a  formulation  does  not  further  simplify  what  has  already  been 
obtained. 


Summary 

The  generalized  forward  map  from  NKOs  to  osculating  modified  equinoctial  elements  described  in  Section  3.1.2 
is  summarized  as: 


m 


VV+Z1) 

A 


f  = 


ZJ2P2^P2  +  z2  -JU 


m2p2^p2  +z2  -p 
A 


cos 

M 

if  z  -  0  a  7=0 

cos 

(q*  -+ 

-  arccos  (^)) 

if  z  cos  j  +  p  sin  j  sin(mt)>0 

cos 

(q*- 

-  arccos  (7)) 

if  z cos  j  +  psin  jsin(mt)  <0 

sin 

(m) 

if  z  =  0  a  7  =  0 

sin 

(q*  + 

arccos  (/)) 

if  z  cos  j  +  p  sin  jsin(mt 

)>0 

sin 

(q*- 

arccos  (/)) 

if  z  cos  j  +  p  sin  7  sin  ( 67/ 

)<0 

h  = 


k  = 


L  - 


Up2  +  z2  —  (p  cos  j-z  sin  jsin(mt))  p  sin  jcos£  +  z(cos  j  cos  ^  sin  (mt)  + sin  ^  cos  (mt)) 
I  sjp2  +  z2  +  (p cos  j  -  zsin  j sin(mt))  ^jz2  cos2  (mt)  +  (p sin  j  +  zcos  jsin(mt))2 

j\lp2  +z2  —  (p  cos  j-z  sin  jsin(mt))  psin  7  sin  ^  +  z  (cos  7  sin  ^  sin  (^)- cos  ^  cos  (^)) 
j  ^p2+z2  +(p  cos  j-z  sin  j  sin  (mt))  ^z2  cos2  (mt)  +  (p  sin  j  +  zcos  jsin(mt ))2 

mt  if  z  -  0  a  7  =  0 

Q"  +  arccos  (/)  if  z  cos  j  +  p  sin  7'sin(^)>0 

Qj  -  arccos  (/)  if  z  cos  j  +  p  sin  j  sin  (mt)  <  0 


(3.95) 


in  which  Q  and  y  are  defined  in  Eq.  (3.40). 
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3.1.3.  Numerical  Test  Cases 

Three  scenarios  have  been  tested  in  order  to  both  verify  the  validity  of  the  mappings  shown  in  Sections  3.1.1  - 
3.1.2  and  study  the  behavior  of  the  osculating  elements.  The  three  test  cases  have  been  chosen  to  be  feasible  mission 
scenarios  in  the  near  future.  For  this  reason,  only  small -displacement  non-Keplerian  orbits  are  taken  into  account, 
which  can  be  generated  with  current  or  near-term  technology.  All  results  of  the  test  cases  have  been  numerically 
verified  by  means  of  MATLAB  functions. 

In  the  following  sub -sections,  all  the  test  cases  taken  into  account  are  described  and  the  mappings  are  shown  and 
analyzed. 


3.I.3.I.  Test  cases 

1.  Displaced  GEO 

Figure  24  shows  the  orbits  (not  to  scale)  of  all  the  test  cases  taken  into  account  related  to  the  GEO  case.  The  value 
of  the  displacement  is  set  to  35  km  for  both  in-plane  and  out-of-plane  displacements  so  that  the  spacecraft  hovers 
above/below  or  inside/outside  the  GEO  station -keeping  box.12 


Figure  24.  Earth-centered  view  (not  to  scale)  of  the  orbits  of  all  the  case  studies  taken  into  account  for  the  GEO  case. 

The  characteristics  of  the  Keplerian  GEO  are  shown  in  Eq.  (2.92).  The  characteristics  of  the  Type  1  NKO  are  as 
follows. 


Type  1  NKO : 


*1 

A 

j\ 

A 


rGEO  Shl(Ayl,) 

=  ^0COS(A;L) 
=  \j El  VGEO 
=  0 
=  0 


(3.96) 


In  this  case,  two  values  of  AZ  are  considered  so  that  both  positive  and  negative  out-of-plane  displaced  NKOs  are 
taken  into  account.  These  values  are  considered  such  that  the  spacecraft  in  the  non-Keplerian  orbit  hovers  above  the 
GEO  station-keeping  box,  which  is  70-km  high.12 
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AA  =  ±  arcsin 


f  \ 

*1 

-  ±  arcsin 

35  km 

V  rGEO  ) 

V  rGEO  ) 

l 

zx  =  ±35  km 
[px  =42,164.1551  km 


(3.97) 


Two  case  studies  are  considered  for  Type  2  NKOs  as  well. 

z2  =  ±35  km 

Pi  ~  VGEO 

Type 2 NKO : 

h  =  0 

4=0 

Two  case  studies  are  considered  for  Type  3  NKOs.  Only  an  in-plane  displacement  is  considered  in  this  case. 


(3.98) 


Type  3  NKO 
(in-plane  displacement) 


z3  =0 

A  =  rGEO  ±  35  kl11 

m3  =PPi  rr-‘ 

h=  o 
£=0 


r3 

' GEO 


(3.99) 


2.  Displaced  global  positioning  system  ( GPS ) 

Figure  25  shows  the  orbits  (not  to  scale)  of  all  the  test  cases  taken  into  account  related  to  the  GPS  case.  The  value 
of  the  displacement  is  set  to  5  km  for  both  in -plane  and  out-of-plane  displacements  so  that  the  spacecraft  hovers 
above/below  or  inside/outside  the  GPS  spacecraft  for  visual  inspection  of  the  same.  Moreover,  a  distance  of  5  km  can 
be  considered  within  the  range  for  proximity  operations.27 


Figure  25.  Earth-centered  view  (not  to  scale)  of  the  orbits  of  all  the  case  studies  taken  into  account  for  the  GPS  case. 
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The  orbit  taken  into  account  for  describing  the  GPS  is  related  to  the  PRN-06  satellite  and  its  orbital  parameters  are 
accurate  with  respect  to  the  Yuma  Almanac  2016*  (week  0896,  GPS  Time  of  Applicability  405504  s).  PRN-06  has 
been  chosen  because  is  the  one  with  the  lowest  value  of  eccentricity  (  e  =  0.22  x  10“3 )  among  the  ones  reported  in  the 
Yuma  Almanac. 


=  0 


Pops  ~  vgps  ~  26,560.9478  km 


GPS  (PRN-06) : 


m, 


r3 

'gps 


GPS 

jGps  =  55.2885  deg 
4,,.  =77.7881  deg 


(3.100) 


The  characteristics  of  the  Type  1  NKO  are  as  follows. 


Type  1  NKO : 


Zi  =  rn « 


,  sin 


.  cos 


(AA) 

(AA) 


7, 


(3.101) 


' GPS 

55.2885  deg 
£  =  77.7881  deg 

In  this  case,  two  values  of  AA  are  considered  so  that  both  positive  and  negative  out-of-plane  displaced  NKOs  are 
taken  into  account.  These  values  are  considered  such  that  the  spacecraft  in  the  non -Keplerian  orbit  hovers  5  km 
above/below  and  outside/inside  the  GPS. 


A  A  =  ±  arcs  in 


> 

Zi 

=  ±  arcs  in 

^5  km 

V  rGPS  J 

V  rGPS  ) 

z  =  ±5  km 

[p  =  26,5609.4733  km 


(3.102) 


Two  case  studies  are  considered  for  Type  2  NKOs. 


Type  2  NKO : 


=  ±5  km 


m 


■-Wr- 


(3.103) 


2  M  A*/  GPS 

j2  =  55.2885  deg 
4  -77.7881  deg 

Two  case  studies  are  considered  for  Type  3  NKOs.  Only  the  in-plane  displacement  is  considered  for  Type  3  NKOs. 

z3=0 


Type  3  NKO 
(in-plane  displacement) 


A  =  >gps  ±5  km 


*"3  AfAV  ' GPS 

j3  =55.2885  deg 
£  =77.7881  deg 


(3.104) 


Data  available  online  at  https  ://celestrak.com/GPS/almanac/Y uma/20 1 6/  [retrieved  25  October  2016]. 


A.  Peloni^|^T^^p,(^^i^proved  for  pub|jc  re|ease.  distribution  unlimited. 


Page  68 


Osculating  Keplerian  Elements  for  Highly  Non-Keplerian  Orbits  -  Final  technical  report 


3.  Displaced  Sun- synchronous  orbit  (SSO) 

Figure  26  shows  the  orbits  (not  to  scale)  of  all  the  test  cases  taken  into  account  related  to  the  SSO  case.  The  value 
of  the  displacement  is  set  to  1  km  for  both  in -plane  and  out-of-plane  displacements  so  that  the  spacecraft  hovers 
above/below  or  inside/outside  the  SSO  spacecraft  for  visual  inspection  of  the  same.  Moreover,  a  distance  of  1  km  can 
be  considered  within  the  range  for  proximity  operations.27 


Figure  26.  Earth-centered  view  (not  to  scale)  of  the  orbits  of  all  the  case  studies  taken  into  account  for  the  SSO  case. 

A  circular  Sun-synchronous  orbit  with  a  radius  rsso  =  1, 000  km  is  considered.  The  value  of  the  inclination  can  be 


easily  computed  considering  the  effect  of  J2  as  j  =  arccos 


2  co, 


’E _  7/2 

2  'SSO 


SSO: 


Zsso  =  0 

Psso  —  rsso  —  km 

^  S SO  —  4p/rsso 

j sso  =  99.4845  deg 

^  ^sso  ~  0 


The  characteristics  of  Type  1  NKOs  are  as  follows. 


Type  1  NKO : 


A 

®i 


J\ 

A 


rsso  sin  (  A/l ) 
^  cos  (A/l) 

99.4845  deg 
0 


(3.105) 


(3.106) 
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In  this  case,  two  values  of  AA  are  considered  so  that  both  positive  and  negative  out-of-plane  displaced  NKOs  are 
taken  into  account.  These  values  are  considered  such  that  the  spacecraft  in  the  non -Keplerian  orbit  hovers  1  km 
above/below  and  outside/inside  the  spacecraft  in  SSO. 


AA  =  ±  arcs  in 


f  \ 


=  ±  arcsin 


r\ kmA 


f  zx  =  ±1  km 

1a  =  rsso 


(3.107) 


Two  case  studies  are  considered  for  Type  2  NKOs. 


Type  2  NKO : 


=  ±1  km 


m 


=4AFr, 


j2  =  99.4845  deg 

4=0 


(3.108) 


Two  case  studies  are  considered  for  Type  3  NKOs.  Only  the  in-plane  displacement  is  considered  for  Type  3  NKOs. 


Type  3  NKO 
(in-plane  displacement) 


P3  =  rsso±  1km 

m3  =  \[AAsSO 

j3  =  99.4845  deg 


[4=0 


(3.109) 


3.I.3.2.  Results  and  Discussion  -  Signatures  of  Non-Keplerian  Orbit  Families 

A  summary  of  the  key  characteristics  of  the  signature  of  the  three  NKO  families  is  given  in  the  following  sub¬ 
sections,  one  for  each  orbital  element  used.  A  detailed  description  of  the  evolution  of  each  osculating  element  for  each 
case  study  is  not  given  here  for  the  sake  of  conciseness. 

Note  that,  in  the  following  sub-sections,  all  the  sinusoidal  curves  are  described  as  cos  (mt  +  cp)  in  which  <p  is  a 
phase  angle  which  is  a  function  of  inclination  j  ,  RAAN  £ ,  and  osculating  element  taken  into  account.  Furthermore, 
f  represents  the  semi-amplitude  of  the  cosine  wave  related  to  the  element  t  . 


1.  Non-Keplerian  orbit  to  osculating  Keplerian  orbital  elements  map 

Semi-major  axis  and  eccentricity  do  not  show  any  special  signature,  being  both  constant  in  time.  The  only  char¬ 
acteristic  that  has  been  noted  is  that,  in  both  cases,  the  in-plane  displaced  orbits  have  the  largest  difference  respect  to 
the  Keplerian  orbit. 


a(t)  =  a0  (3.110) 

e{t)  =  e o  (3.111) 
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The  first  key  characteristic  is  seen  in  the  mapping  of  the  inclination.  Two  trends  can  be  seen: 


in  which 


U(t)  =  i0  ifz  =  0v  j  =  0 

\i (t )  « i  cos [mt  +  cp)  otherwise 


\\i~j\<  0.002  deg  (GPS  case) 
j|/  -  7 1  <  0.008  deg  (SSO  case) 


(3.112) 


(3.113) 


The  right  ascension  of  the  ascending  node  is  characterized  by  three  separate  trends,  as  follows. 

«(*)  =  £  tfz  =  o 

<  fl(t)  =  fl0  +mt  if  j  =  0 

Q(V) «  Qcos(W  +  <^)  otherwise 

in  which 

-  if  z  <  0 

q0  =  2 

3 

—  7i  if  z  >  0 

12 

and 

f|n  - <  0.003  deg  (GPS  case) 

[IQ  - <  0.008  deg  (SSO  case) 


(3.114) 


(3.115) 


(3.116) 


The  argument  of  perigee  is  characterized  by  two  separate  trends,  as  follows. 

[a){t)  =  a>0  if  Keplerian  orbit  v(j  =  0az^0) 
)  =  a\+tnt  otherwise 

in  which 


con 


0  Keplerian  orbit 

<~^  Type  1  (z  <  0)  a  Type  2  (z  >  0),  cox 

3 

—  n  Type  1  (z  >  0)  a  Type  2  (z  <  0) 


J  0  Type  2  a  Type  3  (/?  >  r) 
|  n  Type  1  a  Type  3  (/?  <  r) 


(3.117) 


(3.118) 


A. 
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Lastly,  two  trends  can  be  seen  for  the  true  anomaly,  as  follows. 


in  which 


f  mt  Keplerian  orbit 
ji90  NKO 


J  0  Type  2  a  Type  3  (p>  r) 
|  n  Type  1  a  Type  3  (/?  <  r) 


(3.119) 


(3.120) 


It  is  important  to  note  that  Eq.  (3.120)  shows  that  a  spacecraft  in  a  non-Keplerian  orbit  can  only  be  either  at  the 
perigee  or  the  apogee  of  its  osculating  Keplerian  orbit,  as  it  was  already  discussed  in  the  simplified  case. 


2.  Non-Keplerian  orbit  to  osculating  modified  equinoctial  elements  map 

Semi-latus  rectum  does  not  show  any  special  signature,  being  constant  in  time.  The  only  characteristic  that  has 
been  noted  is  that,  as  for  semi -major  axis  and  eccentricity,  the  in-plane  displaced  orbits  have  the  largest  difference 
respect  to  the  Keplerian  orbit. 

p(t)=p0  (3.121) 


The  first  key  characteristic  is  seen  in  the  mapping  of  the  in-plane  element /.  Two  trends  can  be  seen: 

f(t)  =  0  Keplerian  orbit 

<  w  *  (3.122) 

[f(t)  ~  f  cos(mt  +  (p)  NKO 

However,  the  magnitude  of  oscillation  for  the  out-of-plane  displaced  cases  is  very  small,  as  shown  in  the  following. 

|7<l(r6  ifz^O 

(3.123) 

[/  <  1CT3  ifz  =  0 


The  in-plane  element  g  is  characterized  by  two  separate  trends,  as  the  previous  element. 

(g  (V)  =  0  Keplerian  orbit 

ig(0~Scos(^  +  ^)  NKO  ^  ^ 

However,  the  magnitude  of  oscillation  for  the  out-of-plane  displaced  cases  is  very  small,  with  values  comparable  with 
those  of  the  other  in-plane  element  /  . 


The  out-of-plane  element  h  is  characterized  by  two  separate  trends,  as  follows. 

[h(t)  =  hKep  ifz  =  0 

1  h(t)  «  hcos(mt  +  (p)  ifz^O 


(3.125) 
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The  magnitudes  of  oscillations  respect  to  the  Keplerian  value  are 


(GEO  case) 

\h-hKep\<io-5 

(GPS  case) 

(3.126) 

\h-hKep\<2xl0^ 

(SSO  case) 

The  out-of-plane  element  k  is  characterized  by  two  separate  trends,  as  the  previous  element. 

[k(t)  =  kKep  if  z  =  0 

}  k(t)  «  k  cos  (mt  +  <p)  ifz^O 

The  magnitudes  of  oscillations  respect  to  the  Keplerian  value  are 


(GEO  case) 

\k-kKep\<io~5 

(GPS  case) 

(3.128) 

\k-kKep\<2xlO~< 

(SSO  case) 

Lastly,  one  trend  can  be  seen  for  the  true  longitude,  as  follows. 

L(t)*L(t0)  +  m  (3.129) 

The  absolute  differences  respect  to  the  Keplerian  value  are 


|  L(t)-LKep\*0 

(GEO  case) 

\L(t)-LKep\  <0.001  deg 

(GPS  case) 

(3.130) 

\L(t)-LKep\  <0.01  deg 

(SSO  case) 
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3.1.4.  Non-Keplerian  orbit  to  Osculating  Orbital  Elements  Map:  Summary 
of  the  Signatures 


The  key  characteristics  that  can  provide  a  clear  and  reliable  indication  that  a  spacecraft  is  being  forced  along  a 
non-Keplerian  orbit  are  listed  here. 

The  signatures  of  NKOs  on  the  osculating  Keplerian  orbital  elements  are  the  following. 


KEP: 


i{t)  «  i  cos  (eat  +  (p)  ifz^OAj^O 
Q(/l)  =  Q0  +  mt  ifz^0Aj=0 

<  Q(/l)  «  £lcos(znt  +  (p)  ifz^OAj^O 
co(t)  =  co0  +  mt  if  z  =  0  v  /'  A  0 

i9(f)  =  i90 


(3.131) 


The  signatures  of  NKOs  on  the  osculating  modified  equinoctial  elements  are  the  following. 


MEE: 


f(t)nfcos(mt  +  <p) 
g(t)*  gcos{mt  +  (p) 
h[t )  «  h  cos  (mt  +  (p) 
k{t) «  kcos(mt  +  (p) 


if  z  ^  0 
if  z  *  0 


(3.132) 
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3.2.  Inverse  Map:  Osculating  Orbital  Elements  to  Non-Keplerian  Orbit 

A  generalized  formulation  of  the  inverse  mappings  has  not  being  found.  Nevertheless,  the  formulation  of  the  in¬ 
verse  mappings  shown  in  Section  2.2  for  the  simplified  model  has  been  expanded.  Starting  from  the  Cartesian  position 
and  velocity  vectors  (Eqs.  (3.2)  -  (3.3)),  families  of  NKOs  are  derived  as  a  function  of  the  known  NKO  elements 
given  by  the  simplified  inverse  mappings. 

All  the  families  of  solutions  for  t0  =  0  are  found  by  solving  the  following  system  of  equations  for  the  unknown 
NKO  elements  [ z ,  A  j,  t ] 


r(z,p,w,j,%,t)  = 


Po 

0 


< 


.M 

0 

&0P0 

0 


(3.133) 


in  which  [z0,  p0,m0f  is  the  set  of  NKO  elements  in  the  simplified  case.  Eq.  (3.133)  has  been  solved  via  Mathematica 
function  Reduce  with  the  following  hypotheses. 


{p,p0,m,tu0}>0 
{ z,  z0 ,  f }  e  M 
7'e[  0,tf) 

£  e  [0,  In) 


(3.134) 


Two  main  groups  of  solutions  have  been  found.  These  are  shown  in  Eqs.  (3.135)  -  (3.136). 
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P  =  Po 

UJ  =  ZUq 

'  j  =  0 

£  G  [0,  2 7t) 

}  =  -%!& 


(3.135) 


C  \  n t  =  —  a  k0  ^ 
L 


Po  A  Po  -  Z2o)  +  Poko  +  Zo  (to  -z)K*  Po  A  (C1 

\  / 


B:kg  ^Oa(CvD) 


V 


c2) 


D  \zu  =  m{)P{)  aL  ^ 


-Po  A  (Po  +  4  )  (p0  +  )  *  Z  (Po^  +  ZoK  )  A  (A  V  D2  ) 


a  j  =  2  tan 


^o+0A(CnvC12)  C2:£ 


z  +  z, 


-o  y 


Cu  :  p  =  k0  a  2m0p0t  +  7rk0  (l  +  4 AT, )  =  0 

:0  a  2m0p0t  +  7rk0  (-1  +  4^ )  =  0 


^12  ’  P  ~ 


=  -  —  a  j  =  -2  tan  1 
2 


Po  +K 

z  +  z( 


v  -  ■  -o  J  ^  v^o  1  aoy 

C2l:  p  =  kQ  a  2m0p0t  +  7rk0  (-1  +  4/T, )  =  0  Dn:  p  =  k0  a  2m0p0t  -  7rk0  (-1  +  4^ ) 

C22  :  p  =  -fc0  a  2m0p0t  +  7ik0  (l  +  4^ )  =  0  D^  \  o-  —kn  a  2sj.  oJ  -  7rk„  (l  +  4 K, ) 


A  _  A  /_  \ 

a(C21vC22)  £>,  y  =  ^-a;'  =  2  tan-1  A(DnvZ>12)  D2  :§  =  ~*j  =  2tml  a(D21v£ 

2  v  Po  ^o  y  ^  V  Po  ^o  > 


-u  / 

- u/-  u-  - u  v  -  —  i  /  -  -^21  ’  P  —  ^0  A  _  nK  ^"^1  )  —  ^ 

)|2  :  p  =  — &0  a  2m0p0t  -  7rk0  (l  +  4^  )  =  0  D22  :  p  —  —k0  a  2m0p0t  -  7rk0  (-1  +  4^ )  =  C 

(3.136) 


in  which 


\h=4rf>+4-z 

gZ 


(3.137) 


Note  that  the  terms  Cl2,C22,Dl2,D22  in  Eq.  (3.136)  are  never  verified  because  both  p  and  k0  are  always  positive.  Therefore,  four  separate  families  of  solutions  exist  for 
the  second  group  of  solutions.  For  the  same  reason,  the  constraint  k0  ^  -p0  shown  within  the  solutions  of  the  family  D  does  not  add  anything  and  can  be  safely  removed. 
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The  families  of  solutions  of  the  second  group  (Eq.  (3.136))  are  analyzed  separately  in  the  following  sub¬ 
sections. 


3.2.1.  Family  1 

The  first  family  of  solutions  of  the  second  group  can  be  summarized  by  the  following: 

p  =  K 

zu  =  —  z&qPq /k0 

<  j  =  2tan~'((p0+k0)/(z  +  z0 ))  (3.138) 

£  =  */2 

2m0p0t  +  7tk0  (l  +  4Kt )  =  0 

Note  that  all  the  elements  depend  on  the  value  of  the  out-of-plane  displacement  z  ,  which  can  be  considered 
as  a  free  parameter  of  the  family.  However,  such  element  must  satisfy  the  following  constraints: 

N  *  r0 

<  \z\*\zo\  (3.139) 

Po  (z2  -  z2o )  +  pIK  +  z0(z0-z)k0- p\  *  0 

Note  that  the  third  constraint  listed  in  Eq.  (3.139)  does  not  add  anything  to  the  previous  two  constraints  and, 
therefore,  it  can  be  removed  without  any  loss  of  generality. 

The  solution  related  to  the  time  t ,  shown  at  the  end  of  Eq.  (3.138),  can  be  rewritten  as 

'=-^k(1+“')=^(1+“')=Kf+H  <3140) 

Using  the  true  anomaly  3  =  zut  instead  of  the  time  t ,  it  is  easily  shown  that  the  parameter  K]  can  be  omitted 
without  any  loss  of  generality.  Therefore,  the  first  family  of  solutions  of  the  second  group  can  be  rewritten,  as  a 
function  of  the  out-of-plane  displacement  z  ,  as 

p  =  K 

tu  =  —GT0p0/k0 

<  j  =  2  tan-1  ((/?„+ £0)/(z  +  z0))  (3.141) 

Z  =  2 

3  =  tt/I 


3.2.2.  Family  2 

The  second  family  of  solutions  of  the  second  group  can  be  summarized  by  the  following: 

P  =  k 0 

ST  —  —  STq/Tq  /kQ 

<  ./  =  -2 tarT1  ((p0  +  k(i )/(z  +  70 ))  (3.142) 

£  =  -tt/2 

2 zuQpQt  +  7rk0  (-1  +  4 Kx )  =  0 

Note  that,  also  in  this  case,  all  the  elements  depend  on  the  value  of  the  out-of-plane  displacement  z  ,  which 
can  be  considered  as  a  free  parameter  of  the  family.  However,  this  element  must  satisfy  the  following  constraints 
(in  which  only  the  first  two  constraints  are  shown,  as  discussed  for  the  first  family). 
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fid  ^  rn 

1  I  ,  (3.143) 

LfI  ^  Fol 

As  discussed  for  the  first  family,  the  second  family  of  solutions  of  the  second  group  can  be  rewritten,  as  a  function 
of  the  out-of-plane  displacement  z  ,  as 

P  =  k o 

ex  =  -m0pjk0 

<  j  =  -2Um-'((p0+k0)/(z+z0))  (3.144) 

£  =  -*/  2 
3  -  -7z/2 

Note  that  the  solutions  of  the  second  family  are  characterized  by 

Pi  ~  P\ 

JJJ  2  =  CJ| 

'  h=~h  (3-145) 

6= -6 
92  = 

That  is,  the  solutions  of  first  and  second  families  describe  the  same  osculating  Keplerian  orbits,  as  it  can  be  verified 
by  substituting  Eqs.  (3.141),  (3.144)  into  Eqs.  (3.2)  -  (3.3).  Considering  the  inclination  j  e  ,  the  second 

family  of  solutions  can  be  omitted  without  any  loss  of  generality. 


3.2.3.  Family  3 

The  third  family  of  solutions  of  the  second  group  can  be  summarized  by  the  following. 

P  =  K 

in  =  zu0p0/k0 

<  j  =  2tan~l((z-z0)/(p0+k0))  (3.146) 

£  =  */  2 

2m0p0t  -  7tk0  (-1  +  4 Kx )  =  0 


Note  that,  also  in  this  case,  all  the  elements  depend  on  the  value  of  the  out-of-plane  displacement  z  ,  which 
can  be  considered  as  a  free  parameter  of  the  family.  This,  however,  must  satisfy  the  following  constraints. 

fM*ro 

\  \  '  w  ,  (3.147) 

{{Po  +^o  ){Po+ko)-z{PoZ  +  Zoko)*0 


Note  that  the  second  constraint  of  Eq.  (3.147)  does  not  have  any  solution  in  M  and,  therefore,  it  can  be  safely 
removed  without  any  loss  of  generality. 

As  discussed  for  the  first  family,  the  third  family  of  solutions  of  the  second  group  can  be  rewritten,  as  a  func¬ 
tion  of  the  out-of-plane  displacement  z  ,  as 

P  =  k  0 

iff  =  iff0p0  /k0 

■  y  =  2tan-1((z-z0)/(p0  +  *o))  (3-148) 

£  =  */2 
3  =  -nf  2 
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3.2.4.  Family  4 

The  fourth  family  of  solutions  of  the  second  group  can  be  summarized  by  the  following. 

P  =  k 0 

tn  —  /k0 

<  j  =  2  tan-1  {(—z  +  z0)/(p0  +£„))  (3.149) 

4  =  ~V2 

2m0p0t  -  7rk0  (l  +  AKX )  =  0 

Note  that,  also  in  this  case,  all  the  elements  depend  on  the  value  of  the  out-of-plane  displacement  z  ,  which 
can  be  considered  as  a  free  parameter  of  the  family.  This,  however,  must  satisfy  the  following  constraints  (where 
only  the  first  constraint  is  shown,  as  discussed  for  the  third  family). 

|z|*r0  (3.150) 

As  discussed  for  the  first  family,  the  fourth  family  of  solutions  of  the  second  group  can  be  rewritten,  as  a  function 
of  the  out-of-plane  displacement  z  ,  as 

P  =  K o 

ex  =  m0pjk0 

.  j  =  2 tan-1  ((-z  +  z0  )/(p0  +  k0 ))  (3.151) 

£  =  -*/  2 
3-n/2 

Note  that  the  solutions  of  the  fourth  family  are  characterized  by 

Pa  ~  Pi 

U7  ^  ZU  ^ 

<  u  =  -h  (3-152) 

=  -S3 

That  is,  the  solutions  of  third  and  fourth  families  describe  the  same  osculating  Keplerian  orbits,  as  it  can  be  verified 
by  substituting  Eqs.  (3.148),  (3.151)  into  Eqs.  (3.2)  -  (3.3).  Considering  the  inclination  j  e  \-n,7r) ,  the  fourth 
family  of  solutions  can  be  omitted  without  any  loss  of  generality. 

3.2.5.  Discussion 

Summarizing,  only  two  families  of  solutions  of  the  second  group  exist  and  they  are  described,  as  a  function  of 
the  out-of-plane  displacement  z  ,  as 

Family  1 
p  =  k0 

UJ  =  -ZB 

<  j  =  2  ta 
%  —  ?r/2 
3  -  jz/z 

in  which  j  g  [-7T,  tt)  . 


Family  3 


7oPo  IK 

in'1  ((Po+^o)/(z  +  Zo))’ 


> 

> 


p  =  k0 

tn  =  tu0p0/k0 

7  =  2  tan-1  {{z~z0)/(p0  +k0)) 
%  =  x/2 
3  =  -7r/  2 


(3.153) 
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However,  it  can  also  be  demonstrated  that  the  two  remaining  families  of  solutions  describe,  in  fact,  the  same 
family.  From  Eq.  (3.153),  it  can  be  noted  that  both  orbital  radius  and  RAAN  are  the  same  for  the  two  families. 
Moreover,  both  orbital  angular  velocity  and  true  anomaly  of  the  third  family  are  the  opposite  with  respect  to  those 
of  the  first  family.  In  order  to  have  the  same  non-Keplerian  orbits  from  both  family  1  and  3,  it  must  be 

ri(Zi>7i>') 

_Vl(Wl>0. 

for  any  time  t .  The  result  of  Eq.  (3.154)  is 

jz3=~zi 

(73(23)  =  mod  (7,(2,)  +  n,  22t) 

Because  the  out-of-plane  displacement  z  is  a  free  parameter,  the  first  condition  of  Eq.  (3.155)  can  be  always 
true.  If  the  second  condition  of  Eq.  (3.155)  can  be  demonstrated  to  be  true,  the  equivalence  of  the  two  families  of 
solution  shown  in  Eq.  (3.153)  is  demonstrated  as  well. 

We  wish  to  demonstrate  that 


r3  (Ws.O 
v3(z3,y3,t) 


(3.154) 


7,(2)  j3(~z) 


K 

'2 


From  Eq.  (3.153),  it  can  be  shown  that 

; _ 7,  (2)  _  ( 


a  := 


=  tan 


Po  +^0 

Z+Zg  J 


=  tan 


vx,  +  x2y 


^_AH)=tan-' 

2 


f  ^ 

-z-z0 


Po  +K 


■■  tan 


r  Xx+X2^ 


0 


X 


3  J 


From  the  definition  given  in  Eq.  (3.157),  it  can  again  be  shown  that 


tan  a  =  - 


Xx  +  X2 

O  X{+X2 
tan  B  = - - - - 


tan 


x,  +  x 


-tan 


2  J 


Xx  +  X2 


X 


-a- [5 


3  y 


X, 


tan 


(*-*)-! 


tan  a  -  tan  p  _  Xx+X2  X3 

V  x1  +  x2 


x,  +  x2  x32+(x1  +  x2)2 

x3(x1  +  x2) 


+  tan  a  tan  (3 


1  — 


Ai+A2  X3 


(3.156) 


(3.157) 


(3.158) 


(3.159) 


=>  d-p  =  ^  □  (3.160) 

Eq.  (3.160)  also  demonstrates  that  the  second  condition  of  Eq.  (3.155)  is  true  and,  therefore,  the  two  families 
of  solutions  shown  in  Eq.  (3.153)  are,  in  fact,  the  same  family. 

Because  a  positive  orbital  angular  velocity  is  more  intuitive  than  a  negative  one,  the  third  family  has  been 
chosen  to  fully  describe  the  solutions  of  the  second  group. 

The  two  groups  of  solutions,  for  t0  =  0  ,  are  described  as  follows. 
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P  =  Po 

UJ  —  UJ{) 

7  =  0 

^  e  [0, 2 7r) 


z  =  free  s.t.  \z\  ^  r0 
p  =  k0 

nr  =  m0p0  /kQ 

7  =  2  tan-1  ((z  -  z0  )/(p0  +  k0 )) 

4  =  x/  2 

3  =  -tt/2 


(3.161) 


Note  that  solution  A  shown  in  Eq.  (3.161)  is,  in  fact,  a  particular  case  of  the  more  general  solution  B  .  In 
particular,  zero -inclination  NKOs  are  characterized  by  £  e  [0,2;r)  and  3  =  -<% .  Therefore,  the  general  formula¬ 
tion  of  the  family  of  non-Keplerian  orbits  at  t0  =  0  is  given  by 


z  =  free  s.t.  |z|  ^  r0 

p  =  K 

UJ  =  U70pjk0 

7=2  tan'1  ((z  -  z0  )/(p0  +  k0 )) 
%  =  jr/2 
&  =  -71/2 


(3.162) 


3.2.6.  Generalization  for  any  time 

Equation  (3.162)  shows  the  family  of  solutions  found  for  t0  =  0  .  However,  a  general  expression  for  the  family 
of  NKOs  is  needed  for  any  t0 .  Because  the  NKOs  taken  into  account  are  circular,  it  can  be  chosen  to  use  the 
RAAN  as  a  free  element  instead  of  time  and  leave  the  time  as  t0  =  0  .  The  relation  that  links  the  time  to  the  RAAN 
is 

4=*Vo  (3-163) 

It  has  been  noted  that  the  RAAN  is  the  only  element  that  is  directly  related  to  the  value  of  the  RAAN  of  the 
known  NKO.  Therefore,  the  description  of  the  family  of  NKOs  for  any  time  t0  is  given  by 

z  =  free  s.t.  |z|  ^  r0 

P  =  k 0 

uj  uj0p0  Ao 

(3.164) 

j  =  2  tan'1  [(z- z0)/(p0+k0)j 

£  =  4+^/2 

3  - -nil 
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3.2.7.  Summary 


Given  the  elements  of  an  NKO  characterized  by  the  orbital  plane  parallel  to  the  equator,  the  family  of  all 
possible  circular  NKOs  that  share  the  same  instantaneous  Cartesian  state  with  the  known  NKO  is  given  as  follows. 


z  =  free  s.t.  |z|  ^  r0 

p  =  K 

uj  =  m0p0/k0 

;  =  2  tan'1  ((z  -  z0)/(p0+k0)) 

£  =  4+;r/2 

3  - -7r/2 


(3.165) 


in  which  k0  =  ^pl  +  Zq  -  z  .  In  particular,  zero-inclination  NKOs  are  characterized  by  <Je[0,2 7t)  and 

^  -  4  -  ^  • 


3.2.8.  Numerical  test  cases 


In  order  to  test  the  solution  found  and  shown  in  Eq.  (3.165),  a  numerical  test  case  is  considered.  The  known 
zero-inclination  NKO  is  a  Type  1  NKO  described  by 

Z0  =  rGE0  sin (20 deg) 

Pll  ='G£OCOS(20deg) 

<  ®o  =  \[pI*geo  (3.166) 

Jo  =0 

4=o 

Starting  from  the  particular  NKO  shown  in  Eq.  (3.166),  the  entire  family  of  NKOs  has  been  numerically  in¬ 
vestigated  considering  z  e  [-2 z0 , 2z0  ]  for  two  initial  times  t0 .  For  all  the  test  cases  shown,  the  Cartesian  state  is 
computed,  by  means  of  Eqs.  (3.2)  -  (3.3),  from  the  NKO  elements  found  through  Eq.  (3.166).  In  all  the  test  cases 
shown,  the  error  in  each  element  of  the  Cartesian  position  and  velocity  vectors  is  always  less  than  10"12  km  and 
1CT16  km/s ,  respectively. 
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The  characteristics  of  the  first  test  case  are  shown  in  Eq.  (3.167),  whereas  Figures  27  to  30  show  the  whole 
family  of  NKOs. 


{4=0 

K =° 


(3.167) 


X  [km]  xIO4 

Figure  27.  Test  case  1.  Known  NKO  (black)  and  other  Figure  28.  Test  case  1.  Known  NKO  (black)  and  other 

NKOs,  which  are  part  of  the  family  (red),  are  shown.  NKOs,  which  are  part  of  the  family  (red),  are  shown. 

X-Y  view. 


Figure  29.  Test  case  1.  Known  NKO  (black)  and  other  Figure  30.  Test  case  1.  Known  NKO  (black)  and  other 

NKOs,  which  are  part  of  the  family  (red),  are  shown.  NKOs,  which  are  part  of  the  family  (red),  are  shown. 

X-Z  view.  Y-Z  view. 
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The  characteristics  of  the  second  test  case  are  shown  in  Eq.  (3.168),  whereas  Figures  31-32  show  the  whole 
family  of  NKOs. 


j£o=*/2  ^  J4=° 

\t0  =0  ]tQ  «  6  hours 


(3.168) 


Figure  31.  Test  case  2.  Known  NKO  (black)  and  other  Figure  32.  Test  case  2.  Known  NKO  (black)  and  other 
NKOs,  which  are  part  of  the  family  (red),  are  shown.  NKOs,  which  are  part  of  the  family  (red),  are  shown. 

X-Y  view. 
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3.3.  Impact  of  Noise 

The  impact  of  noise  on  the  NKO  elements  is  now  considered.  This  allows  a  determination  if  the  NKO  signa¬ 
tures  are  still  clear.  The  same  scenarios  used  to  study  the  signatures  have  been  considered.  However,  only  the 
Type  1  NKOs  with  a  positive  displacement  are  used  for  this  part  of  the  study. 

Three  values  for  the  noise  have  been  considered  in  this  study,  consisting  of  upper/lower  boundaries  and  a  value 
inside  this  range.  For  the  upper  boundary,  the  maximum  range  error  |Ar|  =  24  m  and  range  rate  error 

|  Ar|  =  0. 16  m/s  ,  due  to  the  instrumental  errors,  are  considered.28  A  second  set  of  values  for  the  noise  is  considered 

by  dividing  the  upper  boundary  by  a  factor  of  10.  It  will  be  seen  that  these  values  lie  within  the  range  defined  by 
upper  and  lower  boundaries.  A  lower  boundary  on  the  noise  can  be  set  by  considering  the  errors  in  the  orbit 
determination  of  the  test  cases.  Because  these  errors  are  the  result  of  filtering  processes  (i.e.  a  Keplerian  orbit  is 
assumed  in  the  first  place),  these  values  are  considered  here  as  the  lower  boundary,  as  follows.29"33 


(3.169) 


A  Monte  Carlo  analysis  has  been  carried  out  to  obtain  the  maximum  values  of  the  errors  on  the  Keplerian 
inclination  |  A/|  and  RAAN  |AQ|  from  the  aforementioned  errors.  500  samples  of  Cartesian  position  and  velocity 
vectors  in  a  normal  distribution  have  been  considered  with  mean  values  corresponding  to  the  nominal  Cartesian 
states.  The  standard  deviation  considered  for  each  component  of  the  position  vector  is  \\Ar\\ ,  assuming  the  range 

error  being  equal  to  II Aril .  The  standard  deviation  considered  for  each  component  of  the  velocity  vector  is  llAvll , 


assuming  the  range  rate  error  being  equal  to  ||Av|| .  Considering  the  conversion  from  Cartesian  state  to  KEP,  the 


errors  given  by  the  aforementioned  values  provide  an  estimate  of  |  A/|  and  |AQ| .  The  three  sets  of  errors  consid¬ 
ered  here  are,  therefore,  the  following. 


Case  1 

(upper  bound) 


||Ar||  =  24m 
<  |A/|  =  0.01  deg 


180  deg  if  GEO 
0.01  deg  otherwise 


||Ar||  =  2.4  m 

Case  2:  <  | A/|  =  0.001  deg 


if  GEO 


0.001  deg  otherwise 


||Ar||  =  0.1  m 


Case  3 

(lower  bound) 


GEO:  <  |Az|  =  10~3  deg 

|AQ|  =  2x10-3  deg 

iMI  =  lm 

GPS/SSO:  l  | A; |  =  1(T4  deg 
|AQ|  =  KT4  deg 


(3.170) 
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It  is  important  to  note  that  |AQ|  =180  deg  in  the  cases  1-2  related  to  GEO.  This  is  due  to  the  planar  nature  of 
GEO  and,  therefore,  the  fact  that  RAAN  is  arbitrarily  set  to  zero.  Even  a  small  error  can  introduce  a  small  incli¬ 
nation  and,  therefore,  a  value  of  Oe  [0,  2tt)  .  On  the  other  hand,  the  value  of  the  error  on  RAAN  for  case  3  is 
significantly  lower  because  of  the  filtering. 

The  errors  j||Ar||,|A/|,|AQ|}  are  used  to  estimate  the  noise  to  be  considered  for  the  NKO  elements.  A  Monte 

Carlo  analysis  has  been  carried  out  with  500  samples  of  NKO  elements  in  a  normal  distribution.  The  mean  values 
considered  are  the  nominal  values  of  NKO  elements  [z0,p0,  j0,^0} ,  whereas  the  standard  deviations  are 

|||Ar||sin(z0/yO0),^||Ar|,|A/|,|AQ|| ,  respectively.  The  value  of  the  angular  velocity  is  simply  given  by 
m  =  J/u/ p3  ,  in  which  p  is  the  NKO  orbit  radius  with  noise. 

Figures  33  and  34  show  both  GEO  and  Type  1  NKO  in  the  presence  of  noise  as  from  Eq.  (3.170)  for  Cases  1 
and  2,  respectively.  Both  figures  have  been  generated  using  the  Cartesian  position  vectors  resulted  from  the  oscu¬ 
lating  KEP  obtained  through  the  forward  map  from  NKO  elements  with  noise.  The  figures  show  that,  by  looking 
at  their  orbits,  the  Keplerian  GEO  and  the  NKO  are  clearly  distinguishable  for  the  value  of  noise  related  to  Case 
2.  On  the  contrary,  the  same  is  not  true  for  the  three-dimensional  orbits  of  both  the  GPS  and  SSO  cases.  That  is, 
in  the  GPS  and  SSO  cases,  one  cannot  understand  if  the  orbit  they  are  looking  at  is  a  Keplerian  or  a  non-Keplerian 
one. 

Figure  35  shows  a  comparison  between  all  the  scenarios  and  noise  levels  taken  into  account  in  this  study.  The 
evolution  of  the  out-of-plane  MEE  k  has  been  chosen  for  the  comparison.  From  the  plots,  it  can  be  seen  that  the 
first  level  of  noise  is  too  large  to  be  able  to  distinguish  Keplerian  from  non-Keplerian  orbits.  On  the  other  hand, 
with  the  third  level  of  noise,  which  has  been  considered  as  the  lower  bound,  all  the  test  cases  show  a  clear  division 
between  the  two  types  of  orbits.  If  the  second  level  of  noise  is  considered,  the  signatures  of  a  spacecraft  being 
forced  along  a  NKO  are  still  clear  in  the  GPS  and  SSO  cases.  In  fact,  the  constant  trend  of  k  that  is  peculiar  of  a 
Keplerian  orbit  is  clearly  opposed  to  the  sinusoidal  trend  of  k  ,  which  is  characteristic  of  a  NKO.  On  the  contrary, 
the  GEO  case  is  still  too  noisy  to  be  able  to  separate  the  two  trends.  This  is  due  to  the  very  large  error  on  RAAN 
considered  (Eq.  (3.170)). 


- GEO 

Type  1  NKO  -  Out-of-plane  displaced  GEO 


Figure  33.  Case  1.  GEO  and  NKO  with  noise. 


Figure  34.  Case  2.  GEO  and  NKO  with  noise. 


0 

X  [km] 
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Orbit  fraction 


a)  GEO.  Noise  level  case  1 . 


d)  GPS.  Noise  level  case  1. 


g)  SSO.  Noise  level  case  1. 


Orbit  fraction 


b)  GEO.  Noise  level  case  2. 


e)  GPS.  Noise  level  case  2. 


h)  SSO.  Noise  level  case  2. 


c)  GEO.  Noise  level  case  3. 


f)  GPS.  Noise  level  case  3. 


i)  SSO.  Noise  level  case  3. 


Figure  35.  Evolution  of  k  over  one  orbit.  All  mission  scenarios  and  noise  levels  are  shown. 


Lastly,  it  is  interesting  to  study  the  thrust-induced  acceleration  needed  to  follow  the  osculating  elements  in  the 
presence  of  noise.  However,  the  formulation  of  the  acceleration  as  a  function  of  the  osculating  elements  shown  in 
Eqs.  (2.152)  -  (2.157)  has  been  derived  for  a  reference  frame  so  that  viewed  from  an  inertial  frame  of  reference, 
such  solutions  correspond  to  circular  orbits  displaced  above  the  central  body.2  That  is,  the  reference  frame  should 
be  rotated  following  inclination  and  RAAN  of  the  nominal  Keplerian  orbit  before  computing  the  magnitude  of 
acceleration.  Assuming  a  knowledge  of  the  nominal  Keplerian  orbit,  the  reference  frame  can  be  rotated  by  the 
RAAN  and  inclination  following  the  rotations  R3(Q)^R1(/).  This  is  a  good  assumption  since  the  data  are 

noisy  around  the  nominal  value.  Therefore,  the  values  of  RAAN  and  inclination  can  be  considered  as  the  mean 
values  of  the  data  available.  This  way,  the  rotated  Cartesian  state  describes  a  noisy  orbit  with  approximately  zero 
inclination.  Therefore,  the  analytical  formulation  of  the  acceleration  given  either  osculating  KEP  or  MEE  shown 
in  Eqs.  (2.152)  and  (2.155)  can  be  used  to  have  a  first-guess  approximation  of  the  magnitude  of  the  acceleration 
needed  to  follow  the  observed  NKO.  Moreover,  the  value  of  the  magnitude  of  acceleration  can  be  a  further  indi¬ 
cation  of  the  use  of  a  NKO.  Figures  36  and  37  show  the  magnitudes  of  acceleration  computed  considering  the 
second  level  of  noise  for  the  GEO  and  GPS  cases,  respectively. 
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Orbit  fraction  0rbit  fraction 

Figure  36.  Magnitude  of  acceleration  over  one  orbit.  Figure  37.  Magnitude  of  acceleration  over  one  orbit. 
GEO.  Noise  level  case  2.  GPS.  Noise  level  case  2. 


In  Table  10,  the  impact  of  noise  on  the  mappings  is  summarized.  Those  signatures  that,  despite  the  noise,  can 
provide  a  clear  and  reliable  indication  that  a  spacecraft  is  being  forced  along  a  non -Keplerian  orbit  are  highlighted 
for  each  case  study  under  consideration. 


Table  10.  Summary  of  impact  of  noise  on  the  NKO  signatures. 


Signatures 

Issues 

GEO 

GEO  and  NKO  are  clearly  defined  from  the  fol¬ 
lowing  (Cases  2-3  only): 

•  3D  orbits 

•  Osculating  inclination 

•  a  ||  (simplified  analytical  approach) 

•  Case  1  is  too  noisy 

•  From  {Q,  cq,  *9}  and  all  MEE,  GEO  and 

NKO  are  not  clearly  defined 

GPS 

sso 

GPS/SSO  and  NKO  are  clearly  defined  from  the 
following  (Cases  2-3  only): 

•  Elements  {/,Q,/z,k} 

•  a  ||  (simplified  analytical  approach) 

•  Case  1  is  too  noisy 

•  From  {a,#},  GPS/SSO  and  NKO  are  not 
clearly  defined 

•  The  magnitude  of  the  in-plane  MEE  is  too 
small 

•  GPS/SSO  and  NKO  are  not  clearly  defined 
from  the  3D  orbits 

•  ||a ||  can  be  used  only  if  i  and  Q  are  known. 
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4.  Conclusions 

This  work  has  presented  a  mapping  between  highly  non -Keplerian  orbit  (NKO)  geometry  and  classical  orbital 
elements  for  both  the  direct  and  inverse  problem.  Three  sets  of  elements  have  been  discussed  and  mappings  have 
been  derived  in  closed,  analytical  form.  Two  models  have  been  considered  for  the  NKOs:  a  simplified  model,  in 
which  the  NKO  orbital  plane  is  parallel  to  the  equatorial  one,  and  a  generalized  model,  in  which  no  assumptions 
are  made. 

For  the  simplified  model,  it  has  been  shown  that  the  main  drawback  of  using  the  classical  Keplerian  elements 
(KEP)  is  due  to  the  piecewise  formulation  required  for  a  forward  mapping.  A  simpler  formulation  characterizes, 
instead,  a  forward  mapping  that  uses  the  augmented  integrals  of  motion.  However,  the  main  drawback  of  using 
this  set  of  elements  is  the  extreme  sensitivity  to  uncertainties  in  eccentricity  and  true  longitude  that  characterizes 
the  inverse  mapping.  Lastly,  the  modified  equinoctial  elements  (MEE)  provide  both  a  simple  formulation  and 
have  low  sensitivity  to  uncertainties  in  both  forward  and  inverse  mappings. 

For  the  generalized  model,  the  forward  map  has  been  derived  in  closed,  analytical  form.  It  has  been  shown 
that  both  KEP  and  MEE  are  characterized  by  piecewise  formulations.  Three  test  cases  have  been  chosen  that  show 
a  broad  range  of  applicability  of  the  maps.  Key  characteristics  of  the  signature  have  been  highlighted  which  can 
provide  a  simple  and  reliable  indication  that  a  spacecraft  is  being  forced  along  a  non -Keplerian  orbit.  Moreover, 
starting  from  the  Cartesian  state  of  NKOs,  with  their  orbital  plane  parallel  to  the  equatorial  plane,  all  the  families 
of  NKOs  which  have  no  assumption  on  their  orbital  plane  have  been  derived  in  a  closed,  analytical  form  as  a 
function  of  the  out-of-plane  displacement. 

Using  the  same  test  cases  introduced  for  the  study  of  the  generalized  forward  maps,  noise  has  been  added  to 
the  initial  NKO  elements  to  assess  the  impact  on  estimates  of  the  classical  orbital  elements.  It  has  been  demon¬ 
strated  that,  despite  the  noise,  signatures  exist  that  can  provide  a  clear  and  reliable  indication  that  a  spacecraft  is 
being  forced  along  an  NKO.  Lastly,  it  has  been  shown  that  a  first-guess  approximation  of  the  magnitude  of  the 
acceleration  needed  to  follow  the  observed  NKO  can  be  obtained.  Both  KEP  and  MEE  have  equal  advantages  and 
drawbacks  in  the  case  of  the  generalized  forward  map.  Therefore,  in  order  to  clearly  distinguish  a  Keplerian  orbit 
from  a  NKO  using  their  orbital  elements,  it  is  important  to  look  at:  a)  three-dimensional  orbit;  b)  osculating  i  and 
Q  ;  c)  osculating  out-of-plane  MEE  h  and  k  ;  and  d)  magnitude  of  the  thrust-induced  acceleration.  This  can 
provide  an  initial  understanding  for  future,  detailed  analysis  of  optimal  estimation  using  statistical  filtering. 
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